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PREFACE 


This  report  covers  the  analytical  modeling  of  hydraulic  actuator  systems. 

It  also  details  control  design  methodologies.  This  report  may  be  difficult 
to  read  for  the  casual  reader.  It  is  written  assuming  the  reader  has  some 
fundamental  background  in  control  theory. 

I  would  like  to  thank  Contraves  Goerz  Corporation  for  their  notes  descib- 
ing  the  hydraulic  system  which  will  be  used  on  the  Turret  Motion  Base  Sim¬ 
ulator  (TMBS) .  Contraves  Goerz  Corporation  is  the  prime  contractor  for  the 
TMBS. 
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1 , 0  INTRODUCTION 


This  report,  prepared  by  the  Systems  Simulation  and  Technology  Division, 
of  the  U.S.  Army  Tank-Automotive  Command  (TACOM)  describes  an  analytical 
study  of  hydraulic  systems  used  for  laboratory  testing.  Motion  base  sim¬ 
ulators  are  used  by  the  Army  to  test  vehicles  or  vehicle  stibsystems  for 
various  forms  of  performance  and  structural  integrity.  This  form  of 
testing  is  performed  by  hydraulic  systems  which  produce  forces  applied  to 
the  vehicle  simulating  operation  over  terrain  profiles. 

A  considerable  effort  is  being  made  to  develop  analytical  models  of 
laboratory  simulators  so  that  a  complete  assessment  can  be  made  before 
testing  is  conducted.  Analytical  studies  may  point  to  various  problems  with 
the  system  and  control  design  considerations  may  be  made  early  in  the 
development  stage.  The  results  of  selected  command  signals  representing 
terrain  profiles  can  be  simulated  with  the  analytical  model  before  being 
applied  to  the  actual  system.  As  the  level  of  complexity  in  laboratory 
testing  is  increased  the  more  important  becomes  analytical  studies. 

The  Turret  Motion  Base  Simulator  (TMBS)  will  be  a  very  complex  system  used 
for  laboratory  simulation  testing.  The  TMBS  is  a  computer  controlled  6 
degree  of  freedom  motion  base  simulator  which  is  capable  of  driving  an  Ml 
turret.  An  effort  is  being  made  to  model  the  entire  TMBS  system.  The 
starting  grounds  is  to  develop  a  model  of  a  hydraulic  actuator  system  and  a 
means  of  control  design  and  evaluation  of  performance.  Although  the  model 
used  in  this  report  describes  the  hydraulics,  it  only  contains  a  single 
actuator  while  the  TMBS  consist  of  a  six  actuator  system.  Therefore,  the 
results  in  this  study  may  not  describe  the  performance  of  the  TMBS.  At 
this  initial  stage  of  the  TMBS  development,  the  material  presented  in  this 
report  should  only  be  considered  a  study  of  a  general  hydraulic  system.  The 
analytical  study  that  more  directly  relates  to  the  TMBS  will  be  addressed 
at  a  later  time  when  the  entire  system  is  modeled.  A  portion  of  a  model 
describing  the  dynamics  and  kinematics  for  the  TMBS  has  already  been 
established.  The  results  of  this  study  will  be  be  incorporated  in  the 
model  in  the  near  future  so  that  a  model  will  be  constucted  that  contains 
the  six  actuator  systems  coupled  with  the  platform  dynamics. 

This  is  a  first  attempt  in  modeling  hydraulic  systems  in  detail.  An 
earlier  attempt  was  made  in  creating  an  empirical  model  based  strickly  from 
test  data  which  is  described  in  reference  1.  This  technique  consisted  of  a 
curve  fit  model  of  a  hydraulic  system  derived  from  a  measured  closed  loop 
frequency  response  of  the  system.  Although  this  technique  had  resulted  in 
prediction  of  actuator  motion  to  some  degree,  it  was  only  good  for  a  system 
with  no  modifications  or  complexity. 


2.0  OBJECTIVE 


This  report  contains  the  initial  stage  of  an  analytical  study  of  the  TMBS. 
The  primary  objective  is  to  develop  a  hydraulic  actuator  model  along  with  a 
control  design  which  produces  reasonable  performance.  It  is  expected  that 
the  actuator  system  will  have  a  performance  bandwidth  of  10  Hz  with  an 
adequate  stable  transient  response.  The  initial  step  in  achieving  this 
objective  is  to  find  a  means  of  designing  a  control  system  by  using  the 
mathematical  model.  A  methodology  was  developed  which  consisted  of 
writing  a  FORTRAN  program  which  is  used  to  design  control  system 
compensation  by  supplementing  it  with  the  computer  mathematical  model . 
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3.0  CONCLUSIONS 


The  methodology  and  software  developed  to  design  control  compensation  seems 
to  work  very  well.  There  is  definitely  a  correlation  between  the  stability 
margins  obtained  in  the  frequency  domain  and  the  transient  response  in  the 
time  domain.  The  model  shows  that  the  actuator  system  can  perform  a  10  Hz 
bandwidth.  The  step  response  simulated  has  very  good  results  in  terms  of 
overshoot  and  settling  time.  It  is  concluded  that  a  control  design  can  be 
obtained  analytically  without  making  a  linear  reduction  of  the  system. 

The  control  design  obtained  through  this  analysis  may  or  may  not  be 
adequate  for  the  actual  system.  This  will  depend  solely  on  how  accurately 
the  hydraulic  model  compares  with  the  real  system.  The  most  serious 
difficulty  encountered  while  modeling  a  hydraulic  system  is  that  several  of 
the  parameters  required  to  describe  the  system  are  unknown.  What  is  more 
important  at  this  time  is  that  the  methodology  is  now  available  for  control 
design  if  an  adequate  model  can  be  obtained.  The  ground  work  established 
in  this  study  can  lead  to  the  analysis  of  other  hydraulic  systems  used  for 
laboratory  testing  or  for  any  other  application.  Deriving  the  frequency 
response  from  a  non-linear  mathematical  model  can  assist  in  control  design 
or  improvement ing  existing  control  designs. 

The  study  presented  in  this  report  is  a  starting  point  for  more  extensive 
analytical  studies  of  hydraulic  systems.  A  model  of  the  TMBS  has  been 
established  containing  the  dynamics  and  kinematics  which  consist  of  the 
configuration  of  the  6  actuator  orientation  and  platform.  However,  this 
model  required  a  working  model  of  a  hydralic  actuator  with  control.  Now 
that  a  working  hydraulic  actuator  model  has  been  established,  the  next  step 
is  to  incorporate  the  hydraulic  model  into  the  TMBS  model.  The  study  will 
lead  to  investigating  the  control  of  the  platform  by  the  six  actuators. 
Various  methodologies  have  been  developed  by  Contraves  Goerz  Corporation 
which  uncouples  the  effects  of  the  six  actuators  driving  the  platform. 

These  methodologies  consist  of  a  form  of  adaptive  control  which  involve 
considerations  of  the  inertias,  mass  and  actuator/platform  orientations  in 
3-dimensional  space.  This  form  of  adaptive  control,  supplements  with  the 
single  actuator  control  to  complete  a  control  design  for  the  entire  system. 
Future  goals  are  to  investigate  these  methods  thoroughly  by  using 
analytical  models  which  may  lead  to  problem  solving  or  possible 
improvements  with  the  TMBS  system. 


4 . 0  RECOMMENDATIONS 


It  is  impdrtant  to  validate  the  hydraulic  model  developed  in  this  study. 

At  this  time  it  can  be  stated  that  the  model  exhibits  behavior  very  similar 
to  other  known  hydraulic  systems.  The  pressures,  flow  rates  and 
actuator/valve  characteristics  of  the  model  are  reasonable  for  this  type  of 
system.  Nevertheless,  there  is  a  degree  of  uncertainty  to  the  parameters 
used  for  the  model.  Contraves  Goerz  Corporation  is  scheduled  to  perform  a 
single  actuator  test  later  this  year.  At  this  time  the  model  can  be 
validated  and  adjusted  if  necessary.  The  analysis  and  control  design  were 
conducted  using  a  small  mass  for  a  load.  This  was  done  so  that  comparisons 
can  be  made  with  the  single  actuator  testing  using  little  to  no  load  as  a 
starting  point. 

The  search  for  an  optimal  control  design  was  not  real  extensive  at  this 
time.  Once  the  model  shows  some  validity  then  an  extensive  optimal  control 
design  can  be  studied.  This  is  only  when  the  unknown  parameters  can  be 
determined  and  the  proper  configuration  is  known.  (There  are  still 
uncertainties  on  the  inner  loop  configuration) .  Once  this  is  accomplished 
a  further  study  can  be  made  on  optimizing  the  controller  for  a  loaded  case 
consisting  of  a  larger  mass.  10 


5.0  DISCUSSION 


The  starting  point  for  modeling  hydraulics  is  to  obtain  a  set  of  equations 
which  best  describes  the  system.  These  equations  contain  many  parameters 
which  reflect  the  characteristics  of  the  system.  These  parameters  are  not 
always  known  and  must  be  estimated  or  guessed  at  some  way  until 
measurements  can  be  made  on  the  system.  It  may  not  be  possible  to  measure 
some  of  these  parameters  directly.  Some  parameters  may  require 
determination  by  a  combination  of  measurements  and  other  known  mathematical 
relationships.  Once  a  set  of  equations  and  parameters  are  obtained,  a 
model  can  be  developed.  There  are  many  forms  of  computer  software 
available  to  simulate  the  model.  The  one  used  for  this  analysis  is  the 
;  Advanced  Continuous  Simulation  Lanquage  (ACSL)  which  basically  simulates 
equations  in  the  time  domain. 

To  achieve  the  results  of  a  system  with  feedback  control,  a  control  design 
must  be  known  or  derived  from  the  model.  A  control  design  is  part  of  the 
system  that  is  modeled  and  comprises  the  loop  configuration  and  control 
compensations.  For  our  case,  the  control  design  received  was  from  early 
development  and  is  not  considered  acceptable.  A  good  part  of  this  study  is 
to  find  a  means  of  deriving  a  control  design  from  an  analytical  model. 

To  derive  a  control  design,  the  proper  method  using  classical  control 
theory  implies  that  various  frequency  responses  must  be  determined  where  an 
open  loop  response  can  be  evaluated  and  control  compensation  can  be 
designed.  Most  fundamental  courses  in  control  theory  stress  the  importance 
of  linearizing  the  system  equations  (Plant  equations) .  This  is  done  so 
that  the  frequency  response  can  be  easily  obtained  because  the  equations 
reduce  to  transfer  functions  in  simple  Laplace  Transform  format.  With  a 
linearized  set  of  equations  it  is  also  convenient  to  determine  the  eigen¬ 
values  (root  locus  analysis) .  It  is  usually  questionable  when  linearizing 
a  system  if  the  end  results  are  still  truly  representative  of  the  system. 

An  effort  was  made  in  this  study  to  develop  a  software  technique  which  will 
determine  a  frequency  response  from  a  non-linear  mathematical  model 
directly  without  linearizing.  The  primary  philosophy  was  that  no 
questionable  linearizing  techniques  have  to  be  applied  and  the  procedure  of 
algebraically  reducing  the  system  equations  to  a  transfer  function  was  not 
required.  This  method  takes  advantage  of  niomber  crunching  computer  power 
and  is  similar  to  a  frequency  response  technique  which  is  measured  in  a 
laboratory  utilizing  a  signal  analyzer.  The  hydraulic  model  is  then 
evaluated  with  this  technique  and  control  compensation  is  designed  for  each 
feedback  loop.  The  rest  of  this  report  describes  the  details  of  the 
hydraulic  model  and  the  results  of  the  analysis  used  for  this  technique. 


U  HYDRAULIC  MODEL 


Shown  in  Figure  5-1  is  a  diagram  illustrating  the  internal  workings  of  a 
hydraulic  spool /actuator  system.  The  spool  basically  controls  the  fluid 
flow  to  the  actuator  and  is  considered  the  input  state  to  these  equations. 
As  the  spool  is  displaced  up,  the  supply  pressure  creates  flow  to  the 
actuator  bottom  chamber  which  creates  a  pressure  (Pu)  which  in  turn 
produces  an  upward  force  on  the  actuator.  The  negative  actuator  motion  is 
produced  the  same  way  with  a  spool  displacement  down.  The  equations 
describing  the  spool /actuator  system  are  shown  in  Table  5-1.  The  model  is 
a  simplified  case  where  no  overlaps  are  considered  for  the  spool  or  the 
actuator . 
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HYDRAULIC  ACTUATOR  SYSTEM 


SPOOL  ACTUATOR 

A  Xy  f  L,  dUdt 


INPUT:  Xv  Spool  Position 


ACTUATOR:  L  Actuator  Position 

3L/3t  Actuator  Velocity 

Au  ,  Ad  Effective  Actuator  Area 

Vu  ,  Vd  Entrained  Volume 

PRESSURES:  PS  Supply  Pressure 

Pr  Return  Pressure 
Pu  Actuator  Pressure  Up 

Pd  Actuator  Pressure  Down 


FLOWS: 

Qa 

Flow  In  from  Supply 

(Actuator  Up) 

Qb 

Flow  Out  to  Return 

(Actuator  Up) 

Qc 

Flow  In  from  Supply 

(Actuator  Down) 

Qd 

Flow  Out  to  Return 

FIGURE  5-1 

(Actuator  Down) 
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HYDRAULIC  EQUATIONS 

Assume  ideal  Valve  with  Zero  Overlap 


Flow  Area  -  tcDXv  .  W  Xv 


Flow  Q  =  C  Xv  VaP 

Where  C  is  hydraulic  constant 
C  =  Cd  W  V  SQRO 

Ye.iyfi.„..P.ispiaggmgJit 
Valve  Up  Xvu 

Xvu  *  Xv  if  Xv  >  0 
=  0  if  Xv<  0 


D  Spool  Diameter 

Xv  Spool  displacement  from  null 

Q  Flow 

AP  Pressure  across  valve/actuator 
W  Valve  Spool  Circumference 

SQRO  Fluid  Constant 
Cd  Hydraulic  Constant 

Valve  Down  Xvd 

Xvd  =lXvlif  Xv<0 
=  0  if  Xv>0 


Flow  Equations  for  Valve/Spool/Actuator  System 


=  C  Xvu  V  I  Ps  -  Pu  I  Sign  (Ps  -  Pu) 

Qb  =  C  Xvd  V  I  Pu  -  Pr  1  Sign  (Pu  -  Pr) 

Qc  *  C  Xvu  V  I  Ps  -  Pd  I  Sign  (Ps  -  Pd) 

Qd  =  C  Xvd  V  I  Pd  -  Pr  I  Sign  (Pd  -  Pr) 

Actuator  Pressure 

a(Pu)/at  =  B/Vu  (  Qa  -  Qb  -  Au  auat 
a(Pd)/at  =  B  /  vd  (  Qc  -  Qd  -  Ad  ai/at 
where 

Vu*  VuO+  AuL 
Vd  =  VdO  -  Ad  L 


Pu  Actuator  Pressure  Up 
Pd  Actuator  Pressure  Down 
Ps  Supply  Pressure 
Pr  Return  Pressure 
Qa,  Qc  Flow  In  from  Supply 
Qb,  Qd  Flow  Out  to  Return 
(Actuator  Up  Down  Respectively) 

B  Bulk  Moduks  of  Fluid 

^  Au  ,  Ad  Actuator  Piston  Area 
)  (for  up  and  down  respectively) 

VuO  ,  VdO  Initial  Entrained 
Vu  ,  Vd  Entrained  Volume 
(for  up  and  down  respectively) 


Actuator  Force 

F  =  PuAu  -  Pd  Ad 
- - TABLE  5-1 


L  Actuator  Displacement 
aL/at  Actuator  Rate 
F  Actuator  Force 
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The  direction  of  the  valve  displacement  is  detected  and  separated  by  the 
variables  Xvu  &  Xvd  for  up  and  down  motion  respectively. 

This  was  done  so  that  the  up  and  down  cases  can  be  handled  separately  since 
the  two  motions  have  different  characteristics  in  terms  of  flows,  effective 
areas  and  entrained  volumes.  Note  the  u  &  d  subscripts  for  the  remaining 
equations.  There  are  four  relative  pressure  drops  between  the  supply  & 
return  (Ps  &  Pr)  and  the  actuator  pressures  (Pu  &  Pd) .  The  four  relative 
pressure  drops  produce  four  flows  which  are  described  in  the  equations  as 
Qa,  Qb,  Qc  and  Qd.  (See  Figure  5-1)  The  flow  is  assumed  to  be  proportional 
to  the  spool  displacement  (Xv)  times  the  square  root  of  the  corresponding 
pressure  drop  by  a  constant  C.  The  equations  talce  the  absolute  value  of 
the  pressure  drops  since  complex  numbers  do  not  apply  (square  root  of  a 
negative  value) .  The  sign  of  the  pressure  drop  is  considered  by  the  'SIGN' 
function  instead.  Thus  a  negative  relative  pressure  drop  creates  a 
negative  flow.  The  spring  characteristic  of  fluid  is  affected  by  the 
entrained  volumes  (Vu,  Vd,  VuO  and  VdO)  and  the  bul)c  modulus  (fi)  which  is 
assumed  to  be  constant  (no  temperature  variation) .  The  derivative  of 
pressure  is  dependent  on  many  quantities  including  the  actuator  motion.  In 
a  sense  this  can  be  considered  a  natural  internal  feedback  in  the  system. 
Reducing  these  equations  to  a  single  transfer  function  is  very  difficult 
and  is  not  required  for  the  methodology  presented  in  this  analysis.  The 
valve  transfer  function  is  shown  in  Table  5-2  which  describes  the  valve 
dynamics  by  a  second  order.  The  S  represents  Laplace  Transform  notation 
throughout  this  report.  This  equation  describes  the  response  for  a  valve 
where  the  input  would  be  current  and  the  output  is  the  displacement  of  the 
spool.  Also  shown  in  Table  5-2  are  the  various  parameters  used  for  the 
model.  The  most  questionable  and  difficult  to  measure  value  is  the 
constant  C  which  is  the  hydraulic  flow  constant  mentioned  earlier.  The 
various  equations  were  obtained  partially  from  a  combination  of  reference 
2,  3,  and  4.  The  parameters  used  for  this  analysis  are  from  reference  4. 
The  key  element  of  control  is  the  feedback  loops  which  can  now  be 
incorporated  with  these  equations  describing  the  system  to  complete  the 
model.  The  model  was  computer  simulated  using  Advance  Continuous 
Simulation  Language  (ACSL)  which  is  listed  in  Appendix  A. 

It  should  be  mentioned  that  this  hydraulic  model  is  a  simplified 
representation.  There  are  many  nonlinear  hydraulic  fluid  flow  properties 
which  are  assumed  to  be  negligible  in  this  model.  In  addition  this  model 
describes  a  two  stage  servo  system  even  though  it  is  believed  the  TMBS  will 
have  a  three  stage  servo  system.  Since  no  details  have  been  received  on 
the  pilot  valve  characteristics  it  will  be  assumed  that  the  pilot  valve  is 
an  ideal  case.  It  is  also  proposed  that  the  TMBS  will  have  three  valves 
driving  each  actuator  which  are  all  controlled  in  the  same  manner.  This 
was  primarily  done  for  safety  and  backup  considerations.  It  is  also 
assumed  that  this  will  not  change  the  actuator  performance  predicted  by 
this  model.  Some  of  the  assumptions  made  to  simplify  this  model  may  not  be 
valid.  These  modifications  may  be  incorporated  in  the  model  at  a  later 
time. 


5.2  METHOD  FOR  DETERMINING  A  FREQUENCY  RESPONSE 

The  basis  of  designing  a  control  system  using  classical  control  theory  is 
to  obtain  a  frequency  response  of  the  system  (Plant  Response) .  Once  a 
frequency  reponse  is  established,  control  compensation  can  be  designed  to 
desired  stability  margins  and  gain  crossover  (bandwidth)  which  will  reflect 
the  system's  performance  in  the  time  domain.  As  mentioned  earlier,  an 
effort  was  made  not  to  complicate  the  analysis  with  questionable 
linearizing  techniques  which  simplify  but  may  burden  the  accuracy  of  the 
analysis. 
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System  Parameters  Used  In  Modelin( 


Valve  Transfer  Function 


Spool  Position 


Input  Current 


1.5  E-4 


2  mA 

2.533  E-6  S  +  .002228  S  +  1 


Circumference  W  =  3.78  IN 


Actuator 


Effective  Area  =  38.5  IN 


Initial  Entranined  Volume 

3  3 

Up  1697  IN  Down  1336  IN 


Hydraulic  Fluid  Constants 


SQRO  =  1. 


Bulk  Modulus  {  B  )  =  100,000. 


Cd  =  100. 


Pressures 


Supply  =  3000.  PSI 


Return  =  100.  PSI 


Mass 


M  =  --  (Ring  Mass)  =  32.37  Slugs 

6 


Note;  All  frictions  and  hydraulic  leakages  are  assumed  to  be  negligible 


TABLE  5-2 
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The  method  of  determining  a  frequency  response  from  a  non-linear 
mathematical  model  is  shown  in  Figure  5-2.  The  method  is  based  on  taking 
the  Fourier  series  of  both  the  input  and  output  time  histories.  This 
method  can  be  interpreted  as  having  the  analysis  force  the  model  to  behave 
much  like  a  real  system  when  measured  by  a  Fourier  analyzer  in  the 
laboratory.  However  the  input  signal  for  this  case  consists  of  a  sum  of 
sine  waves  at  selected  frequencies  which  span  a  desired  range.  The  Fourier 
coefficients  are  then  used  to  determine  the  appropriate  gain  and  phase  of 
the  response.  The  end  result  of  the  analysis  is  only  considered  at  the 
frequencies  generated  by  the  input  signal.  (A  round  off  technnique  is 
applied  to  round  off  to  the  nearest  milli-Hz)  The  technique  used  in  the 
analysis  was  derived  from  experimentation  with  a  known  linear  model 
describing  a  transfer  function.  This  was  done  so  that  the  correct  results 
can  be  determined  and  compared  with  the  results  of  this  analysis.  Various 
forms  of  input  signals  were  tested  using  this  technique.  An  impulse 
function  was  used  in  several  cases.  Although  the  results  were  close  for 
gain  there  was  an  offset  in  the  frequency  domain  which  was  affected  by  how 
the  impulse  function  was  approximated.  The  phase  could  not  be  determined 
properly  using  this  technique.  Random  noise  was  also  tested  much  like  a 
signal  analyzer  applies  the  input.  This  form  of  input  did  produce  accurate 
results  but  required  a  considerable  amount  of  simulation  time  and  data 
points  to  average  to  a  smooth  spectrum  of  frequencies  (White  Noise) .  The 
sum  of  sine  waves  was  selected  because  it  was  a  much  more  controlled  input 
signal  and  gave  accurate  results.  There  were  only  slight  inaccuracies  at 
the  low  frequencies  (first  few  points)  which  should  not  affect  control 
analysis.  One  drawback  with  this  technique  is  the  limitation  of  data 
points  generated  in  the  frequency  domain.  This  technique  was  used 
previously  to  validate  a  non-linear  math  model  describing  the  Ml  gun/turret 
tracking  system.  The  results  gave  responses  much  closer  to  the  measured 
test  data  than  the  linearized  model  at  various  operating  points. 

Figure  5-3  shows  the  input  signal  used  for  this  analysis  which  consist  of  a 
sum  of  sine  waves  at  the  frequencies  specified.  For  this  particular 
example  the  sine  waves  have  a  magnitude  (or  amplitude)  of  unity.  When  the 
Fourier  series  is  evaluated  at  the  frequencies  generated  the  results  are  a 
magnitude  of  unity  and  angle  of  90°  as  would  be  expected  for  this  case. 

The  input  in  the  time  domain  aproaches  a  function  similar  to  an  impulse 
train  as  more  frequencies  are  added.  One  consideration  when  using  this 
technique  is  determining  what  magnitude  to  use  on  the  sine  wave  generation. 
This  consideration  is  similar  to  determining  what  operating  point  to  choose 
when  linearizing  a  model.  With  some  insight  and  experimentation  one  can 
select  reasonable  magnitudes  to  use  for  a  particular  part  of  the  system. 
Avoiding  saturation  levels  is  the  primary  concern,  however  the  results  are 
magnitude  dependent  as  would  be  expected  for  a  non-linear  system. 

The  frequency  response  analysis  will  be  conducted  at  selected  parts  of  the 
system  model  to  determine  various  plant  responses  and  closed  loop  responses 
for  each  loop  configuration.  The  sine  wave  generated  input  will  be  applied 
to  the  model  at  the  input  of  desired  functions  while  the  output  is  recorded 
for  future  use.  Once  the  data  file  is  created  from  the  model,  a  frequency 
response  analysis  can  be  conducted  using  a  FORTRAN  program  called  Bode. 

This  program  is  listed  in  Appendix  B  and  has  many  applications.  When  the 
frequency  response  represents  a  plant  response  the  program  has  the  feature 
of  designing  compensation  to  the  results.  Since  compensation  is  considered 
a  cascade  function  with  the  plant  response,  the  total  open  loop  can  be 
determined  by  simply  adding  the  gain  and  phase  of  the  functions.  The 
program  has  the  feature  of  doing  this  and  observing  stability  margins. 

The  program  is  menu  driven  and  only  takes  seconds  to  vary  the  compensation 
to  obtain  a  new  open  loop  response.  The  following  section  will  describe 
how  this  technique  was  applied  to  the  hydraulic  actuator  model. 
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INPUT 


OUTPUT 


Fi(t)=  I  Sin(27ift) 


ei{f)  =  TAn\  BNi(f)  /  ANi{f)  )  eo(f)  =  TAN  \  BNo(f)  /  ANo(f)  ) 


PHASE(f)  =  eo(f)  -  0i(f) 

GAIN{f)  =  20Log  (  MAGo(f)  /  MAGi(f)  ) 

This  method  is  evalauted  at  the  frequencies  used  in  the  Sweep  Generator.  f  discrete  frequency 

t  discrete  time 

figure  5-2 
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aPUI; 


INPUT(t)  =  I  A  Sin  (2  7C  Fn  t) 

Where  :  A  desired  magnitude 

Fn  desired  frequencies 

t  simulated  time 

For  our  case  Fn  consist  of  42  selected  frequencies  as  follows: 

Fn  =  .2  ,  .4  ,.6,  .8, 

1.,  1.4,  1.6,  2.0,  2.4,  3.0,  3.4,  4.0,  4.4,  5.0,  5.4,  6.0,  6.4, 
7.0,  7.4,  8.0,  8.4,  9.0,  10.0, 

14.0,  20.0,  24.0,  30.0,  34.0,  40.0,  44.0,  50.0,  54.0,  60.0 
64.0,  70.0,  75.0,  80.0  ,85.0,  90.0,  95.0,  99.0  Hz 

Which  gives  the  following  time  history: 


as. 00 
80.00 
15.00 
10.00 
5.00 
0.00 
-5.00 
-10.00 
-15.00 
-30.00 
-85.00 

0.00  0.50  1.00  1.50  8.00  8.50  3.00  3.50  4.00  4.50  5.00 
_  ()  vs.  TIME  (SECONDS  ) 


SUM  OF  SIN  UAUES  INPUT  MAGS  OF  1 


FIGURE  5-3 


Su3.  ££[liX£QL  DESIGN  MALY  SIS 

Shown  in  Figure  5-4  is  a  block  diagram  of  the  actuator  system.  The  system 
comprises  three  feedback  loops  for  control.  The  states  containing  feedback 
are  position,  rate  and  pressure.  Some  notes  received  show  the  inner  loop 
to  be  force  rather  than  pressure.  This  should  only  result  in  a  difference 
in  gain  in  the  inner  loop.  When  designing  control  compensation  with 
classical  control  theory  most  fundamental  text  books  describe  the  following 
simplified  system: 


Standard  Single  Feedback  Loop 


Where  the  compensation  is  the  only  portion  of  the  system  which  is  subject 
to  change.  The  plant  is  considered  that  portion  of  the  system  which  is 
left  unchanged.  The  design  consist  of  deriving  compensation  which  produces 
the  end  result  or  performance.  For  this  analysis,  the  multi-loop  system 
will  be  reduced  to  three  systems  with  a  single  feedback,  such  as  the  system 
shown  above.  This  procedure  will  be  conducted  one  step  at  a  time  until  the 
entire  system  is  established.  Figure  5-5  shows  the  procedure  in  more 
detail  for  each  step.  The  frequency  response  mentioned  earlier  will  be 
used  as  the  tool  to  evaluate  a  plant  response,  closed  loop  response  and 
design  compensation  by  observing  the  open  loop  response. 


5.3.1  PRESSURE  LOOP  RESPONSE 

Shown  in  Figure  5-6  is  the  inner  loop  block  diagram  which  illustrates  in 
more  detail  the  inner  loop  configuration  including  the  valve  dynamics.  The 
blocks  after  the  compensation  are  considered  the  plant  for  this  case.  The 
plant  response  is  shown  in  Figure  5-7  for  gain  and  phase.  The  valley  at  .4 
Hz  may  be  a  numerical  problem  and  will  not  have  any  influence  on  stability. 
The  peak  at  10  Hz  represents  the  hydraulic  resonance  due  to  compression  of 
the  fluid.  Figure  5-8  shows  the  pressure  compensation  frequency  response 
which  was  derived  by  experimentation  by  observing  the  open  loop  response 
which  resulted  in  Figure  9  for  the  compensation  described.  The  goal  was  to 
establish  the  highest  gain  cross  frequency  with  still  maintaining  adequate 
stability  margins  (See  Figure  5-9) .  The  end  result  is  the  pressure  loop 
response  shown  in  Figure  5-10.  The  closed  loop  response  indicates  a 
bandwidth  greater  than  62  Hz  and  a  very  quick  step  response.  This  gives  an 
indication  to  how  quick  the  system  will  build  up  pressure  to  the  actuator 
for  a  desired  1500  PSI  command.  The  step  response  was  obtained  by  running 
the  computer  model  in  closed  loop  form  and  compensation  design  intact  with 
a  step  command  as  pressure  command.  This  verifies  the  design  which  was 
derived  in  the  frequency  domain. 
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ACTUATOR  SYSTEM  BLOCK  DIAGRAM 


Ol  $ 

O  o; 

O  « 


Actuator  Position 


CONTROL  DESIGN  PROCEDURE 


INNER  LOOP 


SECONDARY  LOOP 


OUTER  LOOP 


FIGURE  5-5 
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PRESSURE  LOOP  BLOCK  DIAGRAM 


Plant  Frequency  Response  of  Pressure  Loop 
Gain 


Phase 


figure  5-7 


Compensation  Frequency  Response  of  Pressure  Loop 


Open  Loop  Frequency  Response  of  Pressure  Loop  « 

f'  Gam  Margin  »  10.2  dB 


Gain 


Gain  Cross  35.3  Hz 


Phase  Margin  =  28.9® 


FIGURE  5-9 


Closed  Loop  Response  of  Pressure  Loop 


Gain  Cross  62. 

^  „  CLOSED  LOOP  RESPONSE  OF  INNER  LOOP  (PRESSURE  LOOP) 

Giain  - - — .  '  j  . . . 

i 

10  _| - - - 

0 


-10 - 

“20  I - 1 — r  T-t  I'l  l - 1 - 1 — 1 — r-r-m - 1 1 — i  i  i  i  i 

1  10  10^ 
dAlN-PRinE-FREOUENCV-RESPONSE  (DB)  vs.  FREQUENCY  (HERTZ) 


DELTA-PRESSURE  (PSD  vs.  TIME  (SECONDS  ) 

FIGURE  5-10 
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The  compensation  for  the  inner  loop  was  determined  to  be  a  second  order 
transfer  function.  A  first  order  transfer  function  resulted  in  a  steady 
state  error.  It  was  decided  to  increase  the  type  number  of  the  system 
(Niomber  of  pure  integrators) .  The  result  is  using  two  PI  (Proportional  + 
Integral)  controls  in  cascade  instead  of  a  single  one. 


5.3.2  RATE  LOOP  RESPONSE 

The  same  procedure  is  now  used  for  the  rate  stage  of  the  system  model  which 
is  illustrated  in  Figure  5-11,  The  plant  is  now  considered  to  be  the 
entire  new  inner  (pressure)  loop  which  includes  the  compensation  which  was 
just  created.  Figure  5-12  shows  the  frequency  plant  response  for  the  rate 
loop  The  compensation  was  a  single  PI  control  shown  in  Figure  5-13.  The 
resulting  open  loop  response  is  shown  in  Figure  5-14.  The  results  are 
shown  in  Figure  5-15  with  the  closed  loop  response  and  the  step  response. 
The  gain  cross  is  over  38  Hz  for  the  closed  loop  response.  The  step 
response  indicates  some  ringing  but  settles  within  .12  seconds.  Perhaps 
the  stablility  margins  could  have  been  improved  for  this  stage.  However, 
at  this  time  the  primary  concern  was  to  establish  an  adequate  response  for 
the  overall  system  which  is  in  terms  of  a  position  command  response. 


5.3.3  POSITION  LOOP  RESPONSE 

With  the  new  pressure  and  rate  loops  having  compensation  design,  the  final 
stage  can  be  evaluated,  completing  the  system.  Figure  5-16  shows  the 
position  loop  configuration.  The  plant  now  includes  the  new  pressure  and 
rate  loops,  which  has  a  response  shown  in  Figure  5-17.  The  compensation  is 
not  shown  because  it  simply  has  a  gain  of  31.6  to  establish  good  stability 
margins.  The  total  open  loop  response  is  shown  in  Figure  5-18.  The  step 
response  is  well  damped  and  has  a  settling  time  of  about  .2  seconds  as 
shown  in  figure  5-18.  The  bandwidth  is  10.1  Hz  as  shown  in  closed  loop 
frequency  response  shown  in  Figure  5-19.  There  is  always  a  compromise 
between  stability  margins  and  bandwidth  performance.  The  bandwidth  can  be 
easily  raised  for  this  system  if  desired.  A  simple  increase  in 
compensation  gain  will  result  in  a  quicker  response  but  produce  more 
overshoot  and  ringing.  The  general  consideration  for  this  application  is 
to  have  a  quick  response  with  minimal  overshoot. 


5.4  RESULTS 

Shown  in  Table  5-3  are  the  results  of  the  compensation  design.  These 
results  are  considered  adequate  at  this  time.  It  is  illustrated  that  there 
is  a  correlation  between  the  stability  margins  and  step  response  char¬ 
acteristics,  Control  compensation  can  be  designed  with  this  technique  and 
may  prove  to  be  of  value  if  the  model  is  accurate  or  is  modified  for 
improving  its  accuracy.  Some  of  the  frequency  responses  derived  here  for 
control  design  are  not  always  possible  to  measure  on  the  actual  system. 

This  stresses  the  importance  of  establishing  a  design  analytically  before 
doing  the  fine  tuning  or  "tweeking”  of  the  controls  in  the  laboratory. 
Analytical  design  gives  insight  into  which  direction  to  go  with  various 
gains  of  the  controller. 
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FIGURE  5-1 1 


Plant  Frequency  Response  of  Rate  Loop 
Gain 


PLANT  REPONSE  OF  SECONDARV  LOOP  (RATE  LOOP) 


Phase 


FIGURE  5-12 
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Open  Loop  Frequency  Response  of  Rate  Loop 


Gain 


Phase 


Gain  Margin  =  5.  dB 
Gain  Cross  14.7  Hz 


figure  5-14 
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Closed  Loop  Response  of  Rate  Loop 


Step  Response 


FIGURE  5-15 
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Position  Loop  Block  Diagram 
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FIGURE  5-1 


FIGURE  5-17 


Open  Loop  Frequency  Response  of  Position  Loop 

Gain  Margin  s  11.7  dB 
Gain  Cross  7.1Hz 


Phase  Margin  *  75® 
Phase  -180®  Cross  27.5  Hz 


FIGURE  5-18 


Closed  Loop  Response  of  Position  Loop 
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TABLE  5-3 


Shown  in  Appendix  C  are  the  results  of  a  simulated  step  response.  This 
Appendix  includes  various  simulated  states  of  the  system  such  as  pressures, 
flows,  entrained  volumes,  spool  displacement  etc.  A  step  response  is 
obtained  by  applying  a  position  step  command  to  the  system.  A  Step 
response  is  used  in  the  analysis  to  verify  control  design  and  give  an 
indication  to  the  performance  of  the  system.  This  type  of  command  signal 
may  not  always  be  practical  to  directly  apply  to  real  actuator  systems  due 
to  harshness.  This  is  especially  true  for  a  position  command  signal.  As 
can  be  seen  by  the  plots  that  the  states  are  driven  to  large  values  the 
first  few  milli-seconds  of  the  simulation.  These  values  point  to 
additional  limits  which  should  be  incorporated  in  the  model.  The  servo 
current  has  a  spike  value  of  2000  milliamps  (2  amps)  which  is  an  excessive 
value . 

Investigations  are  being  made  on  filtering  techniques  which  will  filter 
command  signals  before  they  are  applied  to  the  system.  Previous  laboratory 
testing  consisted  of  detrend  and  filtering  processes  which  modify  command 
signals  as  described  in  some  detail  in  Reference  5.  In  addition  low  pass 
filters  are  used  in  the  laboratory  which  essentially  smooth  command 
signals.  This  would  round  off  command  signals  representing  step  functions 
(Filter  high  frequency  content) .  A  command  signal  which  has  been  filtered 
would  reduce  some  of  the  enormous  overshoots  revealed  by  these  simulations. 

Appendix  D  includes  the  simulated  actuator  response  to  various  levels  of 
step  command.  The  response  changes  as  the  level  of  command  is  varied, 
illustrating  the  non-linear  characteristics  of  the  system.  The  control 
design  procedure  introduced  in  this  report  could  be  conducted  at  various 
signal  levels  to  produce  a  compensation  design  for  a  given  signal  level. 
This  form  of  adaptive  control  could  be  incorporated  in  software  to 
interactively  select  the  appropriate  control  compensation  for  a  given 
signal  level.  With  this  design  a  control  design  could  be  made  which  gives 
good  results  for  all  operating  regions  of  the  system. 

Appendix  E  shows  the  actuator  responses  to  a  variation  of  mass.  These 
results  show  that  as  the  actuator  is  loaded  with  mass  the  transient 
response  is  drastically  changed.  In  fact  the  system  becomes  unstable  as 
the  mass  is  increased  beyond  60  slugs.  This  is  due  to  the  influence  mass 
has  on  the  plant  response  characteristics  of  the  system  which  in  turn  has  a 
direct  impact  on  the  stability  margins.  To  some  extent,  a  change  in  gain 
of  the  rate  loop  should  compensate  (adjust)  to  variation  in  mass.  However, 
this  does  not  completely  improve  the  system  response  as  bandwidth 
performance  is  expected  to  decrease.  These  characteristics  stess  the 
importance  of  adaptive  control  for  the  TMBS  where  variations  in  inertia  are 
expected  to  reflect  on  each  of  the  actuator's  stability  and  performance. 
Additional  concern  is  the  coupling  effects  that  each  actuator  will  induce 
on  each  other.  An  analytical  investigation  will  cover  these  problems  in 
the  near  future. 
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APPENDIX  A 

LIST  OF  HYDRAULIC  ACTUATOR  COMPUTER  MODEL 


PROGRAM  SA  CLOSED  LOOP 

•»  W  —  — 

"  SINGLE  UNCOUPLED  ACTUATOR  SYSTEM  MODEL  for  THBS  " 

n  II 

clnterval  c1nt=0.001 
ARRAY  GFREQISO) 

REAL  DOUT(9dOOi 

INTEGER  INDEX, NTABLE  $”USEO  FOR  SAMPLING  INDEX* 

"  The  following  arrays  ■ 

”  are  used  to  describe  the  transfer  function  coefficient* 

ARRAY  SERVN.SERVDISJ 
"  SERVO  TRANSFER  FUNCTION  COEFFICIENTS* 

CONSTANT  SERVN  -  l.SE-4 

CONSTANT  SERVO  «  2.S33E-6  .  .002228  .  1. 

”  INITIAL  ENTRAINED  VOLUME  and  RING  MASS  * 

CONSTANT  V10sl697.  .  V20«1336.  .  RMASS»194.26 
■  MECHANICAL  PARAMETERS  FOR  VALVE" 

*  BULK  -  OIL  BULK  MODULAS* 

CONSTANT  CD^IOO.  .  SQRO^l.  .  BULK  «  100000. 

*  VALVE  SPOOL  CIRCUMFERENCE* 

CONSTANT  W 1-3.787  ,  N2-3.787 

*  UP  -  LOW  CYLINDER  AREA* 

CONSTANT  AREAl-38.5  .  AREA2-38.5 

*  SUPPLY  8  RETURN  PRESSURE* 

CONSTANT  PS-3000.  .  PR-100. 

INITIAL 

INDEX-0 

CALL  READFREQ(6FREQ. 6LEVEL , NTABLE-) 

END  $”  OF  INITIAL  * 

DERIVATIVE 

«  *********  COMMAND  SIGNAL  ******************” 

*  LEI-  LEG  ERROR* 

"  LCl-  LEG  COMMAND* 

*  LI*  LEG  DISPLACEMENT* 

*  m 

*  GENERATE  SIGNAL  FOR  FREQUENCY  ANALYSIS  " 

"  LCl  POSITION  COMMAND  " 

*  • 

PR0CEDURAL(LC1»GFREQ,QLEVEL,HTABLE.T) 

CALL  GENERATE(LC1.6FREQ.GLEVEL,NTABLE,T) 

ENDS*  OF  PRCEDURAL-GENERATE* 

a  a 

LEI  =  LCl  -  LI 

a  a 

a  *********  POSITION  LOOP  COMPENSATION  *********<> 

*  RCl-  COMPENSATED  LEG  ERROR  SIGNAL* 

RCl  -  KP*(LE1-K)MEGP*INTEG(LE1.0.)) 

a  *********  VELOCITY  LOOP  COMPENSATION  *♦♦*******• 

*  LOEl-  LEG  VELOCITY  COMMAND* 

*  LDl-  LEG  VELOCITY  SIGNAL* 

*  PRCl-  PRESSURE  COMMAND  ???????* 

LDEl  -  RCl  -  LDl 

* 

KT  -  60. 

OMEGT  -  62.8 

PRCl  -  KT*(LDE1+0ME6T*INTEG(LDE1,0.)) 

a  *********  pressure  loop  compensation  *********** 
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ACTUATOR: 


Ln 

LOn 


POSITION 

RATE 


liSj 


**********  CONDITION  THE  SPOOL  POSITION  XV  ********* 
**********  FOR  UPPER  AND  LOWER  CONSIDERATION  ******** 


XVOl  for  spool  down 
XVUl  for  spool  up 


SPOOL  DOWN* 


XVDl*  RSWIXVl  .LT.  0. ,-XVl,0.)$" 

XVUl-  RSW(XV1  .6T.  0..XV1.0.)i”  SPOOL  UP  " 

FLOW  COEFICIENT  Cd 

j  =.C*X*SQRT(DELIA_P) 


where 


=Cd*Wl*SQRT(2/DENS) 


W1»PI*D 


*********  flow  from  pressure  drops  ****************** 


QAl  -  CD*SQRO*Wl*XVUl*SQRT 


ABS 

ABS 

ABS 

ABS 


PS-PU1})*SI6N(1.0,PS-PU1 
PUl-PR  *SI6N(1.0,PU1-PR 
PS-POl  *SI6N(i.0,PS-PD1 
PDl-PR  *SI6N(1*0,PD1-PR 


*********  CYLINDER  VOLUME  VOLUn  B  VOLOn  ************** 
VOLUl  «  VIO  *  L1*AREA1 

VOLDl  =  V20  -  L1*AREA2 

*********  derivative  of  PRESSURE  PDOTUn  (•  PDOTDn*****" 
PDOTUl  «  (BULK/V0LU1)*{QA1  -  OBI  -  AREAl*LDlT 
^  POOTOl  ■  (BULK/V0L01)*(QC1  -  QDl  +  AREAZ*L01J 

*4r******  integrate  for  pressure  PUl  B  PDl  *********** 
PUl  *  INTEGfPOOTUl.O.) 

POl  =  INTEG(PD0T01,0.) 


"  *********  DELTA  PRESSURE  ACROSS  PISTON  ************* 

OELTPl  =  PUl  -  POl 

II  n 

«  *i^^*****  actuator  force  Fn  *************************i 

F1=AREA1*PU1  -  AREA2*PD1  -  FRICT  $•  MUST  LOGIC  THE  FRICtlON* 

m  m 

*  ****  SIMPLIFY  FOR  NOW  ACTUATOR  VELOCITY  AND  POSITION  ******** 

*  ACTUATOR  FORCE  Qn" 

m  m 

MA^SsRMASS/6.  $"  Assume  e  single  unconstalned  mess  fbir  hoW” 
LDl  *  (1./MASS)*INTEG(FI,0.)  $"actuator  velocity  Irtches" 

^  ^  LI  =  !NTEG(LD1,0.) 

end  $''of  derivative* 

****  DISCRETIZE  THE  OUTPUT  INTO  ARRAY  FOR  WRITING  OUTPUT  ******* 

O*************************************************************** 

DISCRETE  SAMPLE 

INTERVAL  DTSMP:::. 00 12207 
STEP=. 0012207 
INDEXsINDEX-fl 
D0UT(INDEX2)=L1 
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“"'•{ISItlKOEXa  .SE.  4099) 

END  $*0F  DERIVATIVE" 

II  It 

terminal^  WRITEER0(=STEP,IH0EX2,00UT) 

END  $•OF  TERMINAL" 

END  l"QF  PROGRAM 

SUBROUTINE  REAOFREQ(eFREQ,SLEVEL,NTABUE) 

1  •  Examle:  freqs.dat’) 

REA0(S,140|FlLEH 

,00 

too  SA0(10:«S0ianil»(IIIABl.E) 

““  If^ffiSwilE)  .IT.  0.)60T0  JOO 
NTABLE-NTABLE+1 
GOTO 

!“  FOaIa?!'  EsioR  OPEHIM  FILE') 

““  Sk?5ir’ 

END 

SUBROUTINE  GEN£RATE(SINPUT.6FREQ.GLEVEL.NTABLE.T) 

*  This  ®’l*put*sl8na*'^to’*a'^fr*quency*’ra5pon5e 

:  ISIlSI*;  tS^'dJKrilSe  function  characteristics. 

*  REAL  GFREQ(50) 

|i:S?I;?99E054.,^ 

K>  '‘'i,J5;?isiS?liT  .  SLE»EL.SIH(J*PI*T*OE«EO(1‘» 

100  CONTINUE 

return 

END 

SUBROUTINE  NRITEERD(STEP.NSAMP.OUTPUT) 

®CTlR:8S8o5fHAME(10).U«IT.»AHE(l6).AUOn.l)OI«YO 

SaRACTERM  y 

CHARACTER*!  COMMA.REPLY 


CHARACTER*2  IOPERATE(20) 
DIMENSION  SMAX(18).SHIN(18} 
REALM  SCALE(l6),6FFSETh6), 
REAL  RMS { 16|^SHEAN ( 16] . OUTPU 


_  ,  _ ],0UTPUT(70  . 

INTE6ER*4  START  SAAP  ELIH.END  SAMP  ELIM 
INTEGER*2  ERD  URIT.HDR  UNIT, I1JATA(T6) 
LOG  I  CALM  TIME 
LOGICAL*!  NEWCHAN, RECHAN 
DATA  ERD  UNIT.HOR  UNIT/10,11/ 

ERD  =  '.ERD* 

HDR  s  '.MDR* 


DATA^ 16, 30000) 


1881 

* 

1888 

1889 

1992 
1920 

1993 


* 

1998 


1766 

* 

201 


how  aany  channels  are  output'} 


NAME  for  channel  ',12) 


WRITE! 5, 1881) 

FORMAT!///.'  ENTER 
READ!5.*)NCHAN 

DO  1998  J-1,NCHAN 
WRITE!5.1888)J 
FORMAT!///,'  ENTER  the  LONG 
READ! 5, 1889) LONG  NAME!J) 

FORMAT!A32] 

MRITE!S,1992]J 

FORMAT!//,'  ENTER  the  SHORT  NAME  for  channel 
REAO(5,1920)SHORT  NAME!J) 

FORMATtAO] 

VRITE!5,1»3]J 
FORMAT!//.'  ENTER  the  UN] 

READ!S, 1920)UNIT_NAME!J) 

OFFSET(J)*0. 

SCALE!J)=1. 

CONTINUE 


*.I2) 


for  channel  ',12) 


00  1766 
DATA! 
CONTINUE 


ISMP=1,NSAMP 

1,ISMP)*0UTPUT!ISHP) 


♦  ' 
♦  ' 
+  ' 
♦  ' 
♦  ' 
♦  ' 


MRITE!5,201) 
FORMAT^//,'  Indicate 


how  new  data  file  Is  to  be  stored.',/. 
Integer  !b1nary]',/, 

/  9 


2  byte  -  , - 

1  e  4  byte  floating  point 
2=8  byte  floating  point 
3=8  byte  complex  !blnary 

4  =  16  byte  complex  !r 

5  =  formatted  floatin 


!  binary)',/, 
ng  point.  Th 


e  format  Is  !Nchannels)E13. 


♦6',/,'$Enter  selection  !0-5):!l  CHOSEN  MOST  COMMONLY)  ') 
READ  !5,*)  KEVNUM 

do  not  let  the  user  choose  complex  numbers 


IF  (KEYNUM  .EQ.  3  .OR.  KEYNUM  .EQ.  4)  THEN 
TYPE*,'  ' 

TYPE*. 'Choose  another  format  besides  complex  numbers.' 
ELSE  IF!KEYNUM  .LT.  0  .OR.  KEYNUM  .6T.  5)  THEN 
TYPE*, '  ' 

TYPE*, 'Selection  out  of  range.' 

ENDIF 

* 

*  open  and  create  files 
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*  begin  writing  header  Inforaatlon 

* 


WRITE(5.26S) 

F0RMAt(//.*$ 


♦/.'  ERb  FORMAT  ASSUMED') 
READ(5,267)  ERD  FILE  0 


Enter  nane  of  the  data  file  to  write  to: 


267  FORMAT (A32) 

*  Create  the  two  file  naaes 
'* 

MRITE(5,4446) 

4446  FORMAT!^  Enter  ERD  title?') 

RE AD (5, 4447) ERD  TITLE 

4447  F0RMAT(A80) 

HOR  FILE  0  »  ERD  FILE  0 

CALU  STRrrRIMIHDR  FILE  0. ERD  FILE  0, LENGTH) 

HDR  FILE  0(LENGTH4l:LEFf6tH-f4T  -  HlJR 
ERD“FILE“0(LENGTH+l:LEMGTH+4j  »  ERD 
0PEH(H0R  UNIT.FILE>H0R  FILE  0. STATUS^' UNKNOWN ' , 
^  ♦  FORM*' FORMATTED', RECL»256)“ 

*  WRITE  OUT  HEADER  DATA 

* 

* 

*  UNKNOWN  KNOWNS 

♦ 

DUMMY  *  'EROFILEVl.OO* 

KEYOPT-0 
NLINES*0 
NBIN*-1 
NBIN— 1 
NBYTE=-1 
COPMA* '  ' 

WRITE(h6r  UNIT,270)  DUMMY 
270  FORMATIAl?) 

WRITE(HDR  UNIT, 280)  ERO.TITLE 
280  FORMAT (A80j 

WRITE! HDR  Jn IT , 290 )  NCHAN , COMMA . NSAMP , COMMA , NL I 


WRITE(HDR  UNIT, 290)  NCHAN.COMMA.NSAMP.COMMA.NLINES.COmA.NBIN. 
6  COMMA .  MBYTE .  COMMA ,  KE  YNUM ,  COMMA ,  STE  P ,  COMMA ,  KE  YOPT 

F0RMAT(6(I7,A),Eli.6.A,!7) 

WRITeIhOR  UNIT.300  scale ( 1),( (COMMA. SCALE! J)),J*2, NCHAN) 
F0RMAT(18TEI3.6,A)) 


F0RMAt(18TEI3.6,A)) 

WRITE(H0R  UNIT,3o6)  OFFSET!!), ((COMMA, OFFSET 
WRITE  H0R~UNIT, 310)  (SHORT  NAME(J), 0*1, NCHAN 
F0RMAT!3lTA8)i 

WRITEIHOR  UNIT, 320)  (LONG  NAME(J), 0*1, NCHAN) 
FORMAT(7(X32)) 

NRITE(HDR  UNIT, 310)  (UNIT  NAME(J), 0*1, NCHAN) 


OFFSET! 0)). 0*2, NCHAN) 
.NCHAN) 


*  write  9+  lines  to  file 

* 


TYPE*,'  ' 

IF(NLINES  .EQ.  0)  GOTO  330 
TYPE*, 'Enter  additional  descriptor  lines' 

DO  330  0=1,NLINES 
REA0(5.280)  LONG 

WRITEIHOR  UNIT, 6560)0, SHORT  NAME(0),RMS(0).SMEAN(0) 

FORMAT!'  CHAN  ' ,I2,3X,A8,3X7'  RMS*  *,E15.3,4X.'  MEAN*  ',E15.3) 
CONTINUE 


A-7 


CLOSE (HDR_UN IT) 

*  write  data  to  file 

IF  (KEYNOM  .EQ.  5)  THEM 

0PEN(ER0  unit, Ft LE»ER0  file  0, form* 'FORMATTED', 

+  STATUS-'^UNKNOWN',  RECr«256T 
ELSE 

pPEN(ERD  UNIT,FILE*ERO  FILE  0. FORM*' UNFORMATTED ' , 
■f  STATUS* 'UNKNOWN',  RECLb256) 

ENDIF 

2  2  byte  Integer 

If  IkEVNUM  .EQ.  0)  THEM 
.  DO  350  J*1.NSANP 
WRITE(ERD  UN1T.ERRM06) 

4-  (liFIX(DATA(L.J)).L*l.NCHAN) 

350  CONTINUE 

h 

*  4  bytb  floatlon  -  binary 

ELSE  IF  (KEYNUM  .EQ.  1)  THEN 
00  370  3«1.NSAMP 

WR1TE(ER0 UNIT,ERR=406) 

*  .  .  TDATAa.<I)>L«l,NCHAN) 

370  CONTIMUE 

e 

*  i  byti  f1 bating  -  binary 

Else  ifikeynum  .eq.  2)  then 

DO  390  J*1.NSANP 
WRITE (ERD  UNIT. ERR»406) 

isio  * 

* 

*  forwattad  output 

ELSE 

DO  410  0*1.NSAMP 
WRITE(ERD  UNIT,40S,ERR*406) 

+  ~  (DATA(L,J),L»1,NCHAN) 

410  .  CONTINUE 


05  F0RMAT(ig(E13.6)) 

05  8rITe!5,407)  J-1  .  . 

D7  FORMAT!/,'  There  were  ',110,'  records  Writtin  but  bbfora  the  f1 

•fie  filled  up.',/,'  Change  NSAMP  In  the  header  flia  accbrdingly. ' ) 

,  ENDIF 

35  CLOSE(ERD_UNIT) 


RETURili 

END 
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APPENDIX  B 

LIST  OF  FORTRAN  PROGRAM  FOR  FREQUENCY  RESPONSE  ANALYSIS 


PROeRAM  BODE 

*************4t***********4r*******************«t******4t**********ik******** 

********  Written  by  A.  L.  HELIHSKI*  US  ARMV  TACOH  AMSTA-RY  ***♦*♦*★ 
********  ******** 

***^********^***********************************^******i******4r********** 
^*************  **************ir******* 

******  FouftlER  ROUTINES  were  developed  by  IEEE  "Programs  ********** 

******  fo^  OlgltSl  signal  processing."  Although  any  Fast  ********** 
******  Fourier  program  which  determines  the  Real  and  ********** 

******  Imaginary  Fourier  coefficients  from  a  time  ********** 

******  history  could  be  used*  ********** 

************************************************************ilr*********** 

*  PURPOSE 

* 

*  this  program  1s  used  to  plot  the  frequency  response  of  dither 

*  the  results  of  a  simulated  model  or  by  directly  putting  In 

*  the  polynomial  transfer  function  coefficients,  the  simulated 

*  model  option  must  have  the  ERO  file  of  the  results.  Simuletlon 

*  must  Involve  using  the  sum  of  sine  waves  as  Input.  The  resulting 

*  file  of  the  output  of  the  transfer  function  time  history  must  be 

*  In  ERD  format.  In  addition  the  sum  of  sine  waves  (Input  of  transfer 

*  Tunctlonl  will  be  re-generated  by  giving  the  same  table  of 

*  frequencies  used  In  the  simulation. 

*  In  addition  the  option  Is  available  to  design  cascade  cOmpehsatloh 

*  If  the  prime  frequency  response  Is  describing  the  plant  response 

*  as  shown  In  the  configuration  below.  The  compensation  frequency  response 

*  can  also  be  plotted  separately  or  combined  with  the  Primary  (Plant)  to 

*  generate  a  total  open  loop  response  for  this  option.  An  option  Is 

*  available  to  determine  crossover  points  representing  stability  margins 


-  ,  PLART 

(PRIMARY)  - 

j  P(S) 

« 

TOTAL  OPEN  LOOP  =  C(S)  k  P(s) 


*  OBJECTIVE 

* 

*  the  iaih  objective  to  developing  this  program  Is  to  have  a  a^mple 

*  means  of  observing  a  frequency  response  with  easy  accest  to  changes  In 

*  the  transfer  function.  (Compensation)  It  Is  also  beneficial  In 

*  designing  compensation  for  a  non-linear  plant  model  by  luRpiamentlng 
"  this  program  with  a  time  domain  simulation  software  package 

*  Tike  ACSL  or  MIMIC. 

* 

**j|r************************************^**i1tgk*************iir*********g^*** 

*********************************A*********^*********^^***iiy*gb********^* 

*****ilt****************4r****************^H^*****************^************ 

*  INPUT  SIGNAL  VARIABLES  ^ 

DIMENSION  SINP(IOOOO)  I  TINE  HISTORY  TO  ISE  BENERATEO 
DIMENSION  AMPI(SO)  I  AMPLITUDE 

DIMENSION  OBliSOl  f  DB  VALUE 

DIMENSION  PHAkltSO)  I  PHASE 

DIMENSION  REALKSOOO)  I  REAL  PART 


a  + 

"  COMMAND—. 


open 
.-X  X 


—  COMPENSATION 
I  C(S) 


B-2 


DIMENSION  AIMAGHSOOO)  I  IMAGINARY  PART 
DIMENSION  GFREQ($0)  t  TABLE  OF  FREQUENCIES  USED 
OUTPUT  SIGNAL  VARIBLES 

DIMENSION  SOUT(IOOOO)  I  TIME  HISTORY  RECOROE 
DIMENSION  AMP0{50)  I  AMPLITUDE 
DIMENSION  OBOISO]  f  OB  VALUE 

DIMENSION  PHASEOISO)  I  PHASE 

DIMENSION  REALOTsOOO)  1  REAL  PART 

DIMENSION  AIMAG0(S006)  I  IMAGINARY  PART 

TIME  t  FREQUENCY  VALUES 

DIMENSION  FREQ„TF(S0).FREQ(5000),SIGTIM( 10000) 
TRANSFER  FUNCTION  VARIABLES 

DIMENSION  OB  TF(50), PHASE  TFISO) 

DATA  HOR  INFO  **t************t**************** 

REAL*4  DATA (10000. 6) 

CHARACTER*8  UNIT  NAME(6), SHORT  NAME(6) 

CHARACTER*32  L0Nff_NAME(6) 

POLYNOMIAL  INPUTS 

REAL  N  NAT(20.20) 

integett  nodr(20).nplys.n  deg 

REAL  POLY_N(40) 

REAL  0  MAT(20.20) 

integeh  ooor(20),dplys.o  deg 

REAL  POLY  D(40) 


I  TIME  HISTORY  RECORDED  FROM  ACSL 
I  AMPLITUDE 
I  OB  VALUE 
I  PHASE 
1  REAL  PART 
I  IMAGINARY  PART 


*********** 


COMPENSATION  POLYS 


REAL  C  N  NAT(20.20) 
INTEGER  Z  NO0R(20).( 
REAL  C_POrV_N(40) 


AT(20.20) 
NO0R(20).C  I 
y-N(40) 

AT(20.20) 
DOdR{20).C  I 
i  0(40)  " 


,NPLYS,C_N_DE6 


REAL  C  D  MAT(20.20) 

INTEGER  Z  0OdR{20).C  DPLYS.C  D  DEG 
REAL  C_POrY_0(40)  “ 


INTEGER*#  NPOINT.NSAMP 
CHARACTER*!  REPLY, NORD 
CHARACTER*8  UNITS 
CHARACTER*32  CHANNEL 
CHARACTER*80  NOTE 
CHARACTER*20  FILETIT 

********************************************************************** 


PI»3. 141592654 

SCAL=180./PI 

EPSILONslO.E-30 


I  Used  for  detection  of  20L06(sina11 ) 


******************************************************************* 
*  INITIALIZE  Set  a11  channels  of  data  to  zero 

A****************************************************************** 

DO  789  Jsl, 10000 
DO  790  K»l,6 

OATA(J.K)sO. 

790  CONTINUE 

789  CONTINUE 
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•  * 


bite^ilhe  thi  users  terminal  type  for  later  use. 


ft 

t 


ilRlTE(5,6) 

TOftMAt(///,’  Enter  the  terminal  Identifier  code:’,/, 


ft  '$Enter  Id: 

REA0(S,*)  ITERM 


1  — >  VT240*,/, 

2  -•>  TAB*  / 

3  — >  TEKTftofilX  40XX’,/, 

4  — >  TEKTRONIX  41XX',/, 
’) 


SS66 

7050 


* 

* 

* 

e 


♦ 
+ 
4.’ 
4* 
♦  ’ 
+  ‘ 

+  * 
+  ♦ 

4< 


MRlTEi5,7050) 

FORMAt?////////////////////////////////, 

!  This  program  evaluates  transfer  fucntloh', 

'  frequency  response.',/,'  It  Is  used  to  plot  frequency', 
Response  and  design  compensation.',///. 

Choose  the  desired  analysis  for  the  PRIMARY  frequency', 
response  (Plant):',///,’  (1)  ACSL  sImulatlon-TIME  HISTORY',/, 
Interpolation  will  be  Included',/, 

(Must  know  ERO  file  generated  from  ACSL  and',/, 
frequency  table  file  used  In  ACSL)',//, 

(2)  Enter  Frcquency/Magnitude/Phase  POINTS  by  hend',/, 
Interpolation  will  be  Included*,//, 

(3)  Enter  transfer  funtlon  In  terms  of  POLYNOMIAL',/, 
COEFFICIENTS',///,'  Enter  1,  2  or  3’) 


REA0(5,*)IMENU1 

Libel  channels  to  this  program's  channel 
configuration  for  later  use. 


NCHAN»6 

LONG  NAME 

LONG'~NAME 

L0NG~NAME 

L0NG“NAME 

L0NG~NAME 

LONGINAME 

SHORT  name  I 
SH0RT”NAMEI 
SHORfNAMEi 
SHORT  NAME I 
SHORTNAMEi 
SHORT“NAMEl 


'GAIN  PRIME  FREQUENCY  RESPONSE' 

'PHASF  PRIME  FREQUENCY  RESPONSE’ 

|='6AIN  EOMPENEATION  FREQUENCY  RESPONSE' 
is'PHASE  COMPENSATIOR  FREQUENCY  RESPONSE’ 

=  '6AIN  TOTAL  CASCADE~FREQUENCY~RESPONSE ' 
=’PHASE_TOTAr_CASCADE  FREQUENCY  RESPONSE' 


* 

* 

* 

5290 

5294 

* 


UNIT  NAME 
UNITHNAME 
UNIOAME 
UNIOAME 
UNIT'^NAME 
UNIT"NAME 


=  'DB’ 
='DEG' 
=  '0B' 
='DE6’ 
=  ’08' 
s'OEG' 


Enter  a  note  (title) 

MRITE(5,5290) 

F0RMAT(*  Enter  a  NOTE  for  this  run  ') 
REAO(S,S294)NOTE 
F0RMAT(A80) 
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*********************************************************************** 

*  OPTION  1  fRlQUEHCi  RESPONSE  FROM  TINE  HISTORY  * 

*  This  option  used  generally  for  a  non-  * 

*  linear  nodel  where  a  time  history  can  * 

*  be  generated.  Simulation  must  be  from  * 

*  sum  of  sine  wave  response  * 


* 

* 

* 

* 

* 

* 

* 


Need: 


ERD  file  of  output  time  history. 
(Created  from  ACSL) 

Table  file  which  contains  the 
frequencies  used, to  generate 
the  Input. 


* 

* 

* 

* 

* 

* 

* 


IFdNENUl  .EQ.  DTHEN 
. READ  INPUT  FILE 

Read  the  date  file  containing  the  tine  history  of  the  output 
signal  of  the  transfer  function.  The  transfer  function  (aodel) 

Must  be  executed  by  a  Input  containing  a  sum  of  sine  waves  for  this 
option. 

CALL  READER0(S16T1H.S0UT.NP0INT. CHANNEL, UNITS) 


-  Deterine  Fourier  Coefficients 

Fourier  subroutine  coMputes  the  fourler  coefficients  of  the  output 
signal  In  terms  of  Real  and  Imamlnary  terms  (Cos  ft  Sin)  for  each 
frequency. 


^  CALL  F0URIER(Sr6TIM,S0UT,NP0INT,REAL0,AIMAG0,FREQ,MM) 

**  READ  FREQUENCIES  AND  GENERATE  THE  INPUT  SIGNAL  (The  same 

*  as  It  was  done  In  ACSL.) 

*  - Read  the  same  frequency  table  as  ACSL 

*  Re-generate  the  same  Input  signal  which  was  used  to  excite  the  model. 

*  (Due  to  memory  limitation  It  was  easier  to  re*generate  the  signal  then 

*  to  create  It  In  the  data  file.) 


* 

* 

* 


771 

* 

* 

* 

* 

* 

* 

* 

** 

* 


CALL  READFREQ(6FREQ.GLEVEL.NTABLE) 

-  Generate  Sun  of  Sin  Waves  of  corresponding  frequencies. 

DO  771  II=1,NP0INT 
T=SIGTIM(Iiy 

CALL  GENERATEISINPQ.GFREQ.GLEVEL.NTABLE.T) 

SINP{II)=SINPO 

CONTINUE 

-  Call  Fourier  for  Input  signal 

Fourier  subroutine  determines  the  fourler  coefficients  of  the  Input 
signal.  Theoretically; 

Real  (Cos  Term)  s  0. 

Imaginary  (Sin  Term)  »  6VEVEL 

CALL  F0URIER(SIGTIN.SINP.NP0INT.REALI.AINA6I,FREQ.MM} 

DETERMINE  AMPLITUDE  AND  PHASE  FOR  ONLY  THE  FREQUENCIES  USED 

INDEX’::! 
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DO  8945  INl.HM 

************************************************************************* 


#*illr*<Ar**** 

iiir'kitititifitit 


********* 


****  First  round  off  frequency  to  the  neerest  bIIH-Hz 
****  Nott*  For  exenple  It  evolds  the  problen  of  4.999  Hz 
****  not  being  recognized  es  5  Hz. 
************************************************************************* 

FREQ(II)»INT(FREQ(I1)*1000.*.5)/1000. 

* 


*  Now  detemlne  If  the  frequency  Is  oa 

*  genereted  Input  signs!  to  Include  In 

***M**********************************i 


one  Included  In  the  origins! 
new  srrsy  of  deta. 


* 
* 

********************M****3»*********************** 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 

* 

it 


it 

* 

* 


*  ' 

* 

* 


-  Deterwlne  If  the  frequency  Is  Included  Mewber  of  roster 

If  the  frequency  Is  the  same  as  one  generated  In  the  Input  signal 
(Included  in  roster-  Freq-Table  member  used  In  Input  generation)  then 
create  In  the  record  (Freq-Hag-Phase)  otherwise  go  to  next  frequency. 

IF(FREQ(II)  .EQ.  6FREQ( INDEX) )THEN 

Create  new  frequency  array 

FREQ_TF(INDEX)=FREQ(II) 

Create  aaq>11tude  array  of  Input  * 

AMPI ( INOEX)«SQRT(REALI ( II )**2*AIHA6I ( I I )**2) 

Create  dB  array. 

Avoid  L0610  of  very  small  numbers 


IF(AMPIj INDEX)  .LT.  EPSILON)THEN 


t  -S80  Db  Is  small,  make  as  limit 


DBl(lNDEX)=-580. 

ELSE 

08I(INDEX)>20.*AL0610(AMPI(INDEX)) 

END  IF 

Create  amplitude  and  dB  array  of  output  signal  * 

AMPO (I NDEX 1 =SQRT ( RE ALO ( I I ) **2+AIMAG0 (I I ) **2 ) 
IF(AMPOpNDEX)  .LT.  EPSILONJTHEN,  _ 


I™ 


DBO(INDEX)s-S80. 

ELSE 

^^^080 ( INDEX ) =20 . MLOG 10 ( AMP0{ INDEX) ) 

Create  Phase  array  of  Input  and  output  signal 

lF(REALI{in  .EQ.  0.)THEN 
IF(AIMAGl{in  .LT.  0.)THEN 
PHASEI(IN0EX)=-90. 

ELSE 

PHASEI(INDEX)=90. 

END  IF 

DETERMINE  WHICH  QUADRANT  THE  INPUT  PHASE  IS  IN 


580  Db  Is  small,  make  as  limit 


ELSEIF^REALI(II) 
PHASE I 


:eali(ii)  .gt.  0.  .and.  aim 

=ABS ( AIMAGI ( I I )/REALI ( I I ) ) 
(IN0EX)=  SCAL*ATAN(R, 


AIMAGI(II)  .GT.  0.)THEN 


lATIOI ) 
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* 

* 

* 

* 


* 

* 

* 


* 

* 


ELSEIF(REALI(II) 
RAT10I«ABS 
PHASEKINDEXj 


ELSEIF(REALI 


(II)  .LT.  0.  .AND.  AIN. 
(AIMA6I(II)/REALI(Ilh 
EX)»180.  -  SCAL*ATAN(ft, 


AINAGKII)  .GT.  0 


(REALI(It)  .LT. 
0I-ABS(AIM6I(I 


_  lATIOn 
0.  .AND.  AlHAGI(lf) 


A6I(II)/REALI(in) 
180.  -I-  SCAL*ATAN{R. 


RATlOr»ABS^ 

PHASEIIIND^X 
ELSEIF(REALI{I  . 

RATI0I=:ABS{AINA6I  ( II  )/REALI  ( 1 1)) 
PHASEI(INDEX)=360.  -  SCAL*ATAN(Ri 


)  .GT.  0. 


ATIOI] 


.AND.,AlMA6I(li) 
ATIOI) 


.LT.  0 


.LT.  0, 


* 

* 

* 

* 


ENDIF 
OUTPUT  PHASE 


IF(REALO(II)  .EQ.  0.)THEN 
1F(AImAG0(II)  .LT.  0.)THEN 
PHASEO{IN0EX)s-90. 

ELSE 

PHASE0(IN0EX)=:90. 

ENDIF 

DETERMINE  WHICH  PHASE  THE  OUTPUT  IS  IN 


ELSEIF(REALO{II 
RATIOO=j 
PHASED 
ELSE IF 
RATI 


i®T‘.<L_:A!!O..AIMAGO(lI)  .GT. 

RATIOO) 
GO(II)  .GT. 


bABS (AIMAG0( I I )/REALO (II)} 

>(INOEX}»  SCAL*ATAN( RATIOO) 

aREALO{fn  .LT.  0.  .AND.  AIMAgO(II)  . 
=:ABS(AiM60(II)/REAL0(II)} 
PHASEO(INOEX}»  180.  -  SCAL*ATAN( RATIOO) 
ELSElF(REALO(tn  .LT.  0.  .AND.  AIMAgO(II)  .LT. 
RATIOO»ABS(AiHAGO(II)/REALO(II)) 
PHASEO(INDK)=  180.  f  SCAL*AfAN(RATIOO) 
ELSEIF(REALO(tn  .GT.  0.  .AND.  AImAgO(II)  .LT. 
RATI00=ABS( AimAGO( 1 I )/REALO( I I ) } 
PHASE0(IND»)«  360.  -  SCAL*ATAN(RATIOO) 

ENDIF 

Create  dB  and  phase  array  of  transfer  function 
08  TF(INDEX):=DB0( INDEX) -OBI (INDEX) 

phAse_tf(in6ex)=phaseo(index)-phaSei(index) 

INDEX=INDEXi-l 


■k 

8945 

* 

* 

* 

* 


ENDIF 
CONTINUE 

WHICH  DELTA  F  INTERPLOATE 
WR1TE|5,,4659)FREQ_TF(  1) 


4659  FORMA 
+ 

♦ 


!W»  DELTA-F(In  Hz)  FOR  INTERPOLATION?*, 
Choose  smaller  or  equal  to  ',F10.4) 


REAO(5,*)STEP 


.)THEH 
.)THEN 
.  )THEN 


0. )THEN 
0. )THEN 
0.)THEN 
0.)THEN 
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* 

* 

* 

* 

* 


*** 

* 

* 


44333 

* 


-  Interpolate  date 

Interpolate  the  frequency  response  data,  aore  points  will  be  generated 
for  a  constant  step  (Delta  Hz).  This  Is  done  for  plotting  purposes  used 
later. 

fENO=FREQTF( INDEX-1) 

NSAMPalNDEX-l 

CALL  IHTEP ( 1 , NSAMP , STEP, FREQ_TF, DB_TF , PHASE_TF , DATA) 

NSAHPsNSAMP-1  I  TRYING  TO  ELIMINATE  ONE  MYSTERIOUS  BULLSHIT  POINT 

FORM  TOTAL  OPEN  LOOP 

Add  primary  and  Compensation 

DO  44333  Jsl, NSAMP 

0ATA(J.5}:=DATA(J,1}4DATA(J,3} 

OATA(J,6)=OATA(J,2)+OATA(J,4) 

CONTINUE 


* 

* 

* 

* 

* 

★ 


3000 


3003 

3002 

* 

4669 


* 

*** 

* 


99333 

* 


ENOIF  1  OPTION  IMENU1=1 

OPTION  2  ENTER  DATA  POINTS 
**************************** 

This  option  Is  available  when  only  a  few  frequency  date  points  are 
known.  (Usually  from  test  deta  or  extracted  from  a  plot.)  This  option 
will  Interpolate  points  In  between  the  given  ones.  Be  Sure  to  Include 
enough  points  to  desribe  the  response.  The  more  drastic  the  changes  the 
more  points  required. 

IFdMENUl  .EQ.  21THEN 
HRITE(5,3000) 

FORMAT!'  Enter  how  many  points  to  enter  by  hand', 

*  '  (This  contains  each  Frequency,  Gain  and  PhlSe  as', 

+'  one  point)') 

REA0(5.*)NPTS 
DO  3002  3»1,NPTS 
HRITE(5,3003jj 

FORMAT!'  POINT  ',12,'  Enter  FREQ(Hr),MAG(dB) ,PHASE(Deg) • ) 
READ(5,*)FREQ_TF(J),DB_TF(J),PHASE_fF(J) 

CONTINUE 

HRITEjS, 4669) FREQ  TF(1) 

0ELTA-F(In  Hz)  FOR  INTERPOLATION?', 

<f  '  Choose  smaller  or  equal  to  *,Fi0.4) 

READ(5,*)STEP 

F_END=FREQ_TF( INDEX) 

CALL  INTEP(1,J, STEP, FREQ_TF,DB_TF,PHASE_TF, DATA) 

FORM  TOTAL  OPEN  LOOP 

DO  99333  J=l, NSAMP 

0ATA(J,5)=DATA(J,1)+DATA{J,3) 

DATA{j,6):=DATA(J,2)->-DATA(J,4) 

CONTINUE 

ENDIF  I  OPTION  2  (POINTS  BY  HAND) 
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*M*******M***********«****«**llr********************************* 

***♦  OPTIION  3  ENTER  PRIME  FREQ  BY  POLYNOMIALS  **♦**!»*♦*****♦**♦ 

*  This  option  gives  the  frequency  response  for  a  transfer  function 

*  Model  described  by  polynoalals  In  the  S  (Laplace)  donaln. 
***************************************************************** 


4689 


5689 


* 

* 

* 

a 

a 

a 

a 


a 

aaa 

a 

a 


IFflNENUl  .EO.  3 j THEN 

CALL  ENTER  NP0LYS{N  MAT, NODR, NPLYS.N  OEG)  t  ASK  NUMERATOR 
CALL  ENTER  DP0LYS(0^r,00DR.DPLYS,D  DEG)  !  ASK  DENOMINATOR 
CALL  MULTIFLY  POLYSTH  0E6.H06r,NPLYS7N  MAT, POLY  N)  I  MULTIPLY  NUM 
CALL  MULTIPLY;pOLYS(D3E6,DODR.DPLYS,DIMAT,POLY;D)  !  MULTIPLY  DEN 

WRITE|S,4689) 

0ELTA-F(In  Hz)  FOR  INTERPOLATION?',/, 

•  {  .1  Generally  Used  )M 
READ(S,*lSTEP 

WRITE(5,5689) 

'  hoo.  Generally  Used)*) 

READ{5,*)F_END 

-  Poly  Response 

Detenelne  the  frequency  response  for  the  Nunerator  and  Denoelnator 
Product  Polynoelat  (Poly  N  K  Poly  D  respectively).  Find  response  for 
frequency  range  STEP  to  F  END  by  STEP.  The  channels  for  the  Prloary 
response  are  I  &  2  for  Gain  and  Phase  respectively. 

CALL  POLY  RESPONSE{N_DEG.POLY  N,D  DEG, POLY  D,STEP, 
F_END,l,2,0ATA,llSAMP) 

FORM  TOTAL  OPEN  LOOP 

By  adding  the  Primary  and  Compensation  response.  (Cascade) 

DO  7333  Jxi  nsAMP 

DATA(J,S)>=DATA(J,1}4-DATA(J,3) 

0ATA|J,6)=.0ATA(J,2)+0ATA(J,4) 

CONTINUE 


7333 

* 

* 

ENDIF  I  OF  OPTION  3 

aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

*  2ND  MENU  * 

****************************************************************** 

* 

WRITE (5 ,900) 

U( . 


901 

900 


♦ 

+ 

♦ 

+ 

+ 

+ 

+ 

♦ 


frequency  response’,/, 
COMPENSATION  by  polynomials' 
OR  EDIT  polynomials',/. 


PLOT  frequencj 
Enter 

CHECK  OR  EDIT  polynoR 

Determine  CROSSOVERS  for  total  response',/. 

Enter  NEW  PRIME  FREQUENCY  RESPONSE',/, 

STOP  /  QUIT’,///, 

Select  1,2, 3, 4, 5  or  6') 

READ(5,*jlMENU2 

****************************************************************** 

*  OPTION  1  (2ND  MENU)  PLOT  FREQUENCY  RESPONSE  * 

aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa 

IF(IMENU2  .EQ.  1)THEN 
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*  *  *  *  «  «  « 


CALL  rr  (>L0r(ITERH.N0TE.NSAMP,STEP.NCHAN.0ATA.L0N6  NAHE, 
+  UNIT  NARE) 

GOTO  901 

ENOIF  t  INENU2  1ST  OPTION 

****************************************************************** 

*  OPTION  2  {2ND  MENU)  ENTER  COMPENSATION  POLYS  * 

****************************************************************** 

* 


Enter  compensation  by  means  of  transfer  function  polynomials 
In  S  (Laplace)  domain,  same  technique  used  above  for  option  3 
of  entering  primary  response. 


1F(IHENU2  .EQ.  2)THEN 
CALL  ENTER  NP0LYS(C 


Call  enter  npolys(C  n  hat.c  nodr.c  nplys,c  n  deg) 

CALL  ENTER~DP0LYS(C  D^AT.C  D0DR,C_DPLYS,CD)DEG} 

CALL  MULTIPLY  P0LYSTC~N  OEG,C_HOOR,C_NPLYS,C_N  MAT.C  POLY  N) 
CALL  MULTIPLY“POLYS{C  D'DEG.C  DODR.C  DPLYS.C  DTIAT.C  POLY  D) 
CALL  POLY  RESRONSE(C  R  BEG.CROLY  M.C  D  DEG.t  FOLY  D7 
STEP ,  F“ERD,3, 47DATA7NSARPT 


***  . FORM  TOTAL  OPEN  LOOP 

*  by  adding  primary  and  compensation  reponse  (Cascade) 

DO  9333  J>1.NSAMP 

OATA(J,5}aDATA(J.l}*DATA(J,3) 

OATA(J,6hOATA{J.2)*DATA(J.4) 

9333  CONTINUE 

* 

GOTO  901 

* 

ENDIF  1  INENU2  2ND  OPTION 

************************************************************************ 

*  OPTION  3  (2N0  MENU)  EDIT/LIST  POLYNOMIALS  * 

* 

IF(IMENU2  .EQ.  3)THEN 


IFdMENUl  .EQ.  3)THEN 
WRITE(5,3335) 
FORMAt(///////////// 


’transfer  FUNCTION’,//, 

♦  ’  (2)  COMPENSATION  TRANSFER  FUNCTION’,///, 

♦  •  ENTER  1  or  2’) 

REA0(5.*)IC0MP 

ELSE 

ICONP=2 

ENDIF 

* 

LIST  PRIMARY  POLYS  * 

* 

llr*******************1^*******4r******************«**************** 

IF(ICOMP  .EQ.  DTHEN 
->•  Liat  Individual  polynomials 
WRITE(5,5550) 

CALL  LIST_P0LY(1,N_MAT,P0LY_N,N0DR,NPLYS,N_DE6) 


B-10 


WRITE(5.5S51) 

CALL  LI ST_POLY(  1 , 0_HAT ,  P0LY_0 , DODR ,  DPLYS ,  D_OEG ) 

LIST  PRODUCT  POLYS,  unless  the  seme  as  above 

IF(NPLYS  .ST.  DTHEN 
WRITe75,4441) 

4.  ***********  //\ 

CALL  LIST  PdLH2.N  HAT. POLY  N.NODR.l.N  DEG) 
ENOIF  -  -  - 

IF(DPLYS  .GT.  DTHEN 
WRITE (5,4442) 

▲  •  *********1  //\ 


CALL  LIST_P6lV(2.D_HAT,P0LY  D.OOOR.l.D  DE6) 
ENOIF 


* 

* 


ELSEIF(1C0NP  .EQ.  2)THEN 

LIST  COMPENSATION  POLYS 


NRITE(S,S530) 

CALL  LIST_POLY(l,C_N_HAT,C  POLY  N. 

+  _  _  -  gjiQjjp  c  UPLYS.C  N_OEG) 

WRITE(5,553D 

CALL  LIST  POLYd.C  D  HAT.C  POLY  D, 

^  ♦  C_DODR,C_BPLYS,C_D_DEG) 

*  LIST  PRODUCT  POLYS,  unless  the  same  as  above 

* 

IF(NPLYS  .GT.  DTHEN 
WRITE{S,5534) 

5534  FORHAt  (////////////////////////////////// , 

♦  *********  COHPENSATION  PRODUCT  NUMERATOR  *******•//) 

CALL  LIST  P0LY(2,C  N  MAT.C  POLY  N,C  NODR.l.C  N  DEG) 

ENDIF  "  ___ 

IfIdPLYS  .GT.  DTHEN 
MRITE(5,5537) 

CALL  LIST  P0LY(2,C  D  HAT.C  POLY  D.C  DODR.l.C  D  DEG) 

ENDIF  !  DFLYS  .GT.“l  "  "  ~ 

itifititititit'kititit'k'kitit'kititlt'kltitititititit'kitit'kltltitit'kitifkit'kiUtifkititit'kif'k'kititit'k'klfk'kitit’kitit’k'kitit'kit 

*' 

ENOIF  I  ICOMP  I  or  2 

* 

*********1^****ifr****1^********4r********************************<fr******** 

**  EDIT  POLYNOMIAL  * 


FORMAT 

I  iekhltlrk 
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WR1TE(5  5555) 

S5SS  formaT(>///////////////////////////////////////. 
*  '  Oo  you  want  to  «ake  a  change?  (Y  or  N)') 

READ(5,77777)REPLY 
FORMATiA) 

IF(REPLY  .EQ.  'Y'  .OR.  REPLY  .EQ.  'y’jTHEN 
NUMERATOR  OR  DENOMINATOR  ************** 


nm 

*■ ' 

* 

** 

'it' 

*. . 


66666 


18 

* 


WRITE (5,66666) 

oy, 

REA0(5.18)N0R0 

FORMAT(A) 


* 

iit. 


IFdCOMP  .EQ.  DTHEN 

IF(NORD  .EQ.  ‘N*  .OR.  NORO  .EQ.  'n*)THEN 

k*******************************************************1 

^CHANGf  PRIMt 

::..—  ')CHANeE  PRIME  NUMERATOR . 

GALL  EDIT  POLY(l,N  MAT,NOOR,NPLYS.N  DE6) 
ELSEIF(NORD  7EQ.  *d’  TOR.  WORD  .EQ.  ’iDTHEN 


* 

* 

* 


.  CHANGE  PRIME  DENOMINATOR  . 

CALL  EDIT  POLY(Z.O  MAT, DOOR, DPLYS,D  DEG) 
ENDIF  -  - 

RECREATE  FREQUENCY  RESPONSE  FOR  PRIME  AND  TOTAL 


3333 


CALL  MULTIPLY  POLYS (N_DEG, NODR, NPLYS.N  MAT, POLY  N) 
CALL  MULTIPLY  POLYS (D  DEG.DODR.DPLYS.D  MAT, POLY  D) 
CALL  POLY  RESFONSEIN  DEG, POLY  N,D  DEG, TOLY  D,STrP, 
F  ERD, 1.2, DATA, NSAMP)  "  “  “ 

DO  3333  J=1,NSAMP 
DATA7j,5)=DATA(0, 1)+DATA( J,3) 
0ATA(J,6}=0ATA(J,2)+DATA(J,4) 

CONTINr 


********************************************************** 


I  ICOMP  MUST  EQ  2 
***********&************************************************ 

*  ^HANGt COMPENSATION  POLY  ft  RESPONSE  * 

************************************************************ 


* 

i. 


* 

* 


- .  CHANGE  COMPENSATION  NUMERATOR . 

’  IF(NORO  .EQ.  *N'  .OR.  NORD  .EQ.  *n')THEN 

CALL  EDIT_POLY(l.C_N_MAT,C_NODR,C_NPLYS,C_H_DEG) 

.  CHANGE  COMPENSATION  DENOMINATOR  . 

ELSEIF(N0RD  .EQ.  'D*  .OR.  NORD  .EQ.  'd'lTHEN 

CALL  EDIT_P0LY(2,C_D_MAT,C_D0DR,C_DPLYS,C  0  DEG) 
ENDIF 

PECREATE  FREPUENCY  RESPONSE  FOR  COMPENSATION  AND  TOTAL 
CALL  MULTIPLY  POLYS (C  N  DEG.C  NODR,C  NPLYS, 

+  "  c_h_mat,i:_poly_R) 
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14  ilWWilyiiMylw  Idy  iitf *4U3I  wi 


IF(1MENU2  .EQ.  4)THEN 
CALL  CROSS.OVER ( 5 , RSAMP , STEP , DATA ) 

GOTO  901 
ENOIF 

>  OPTION  5  (2N0  MENU)  RETURN  TO  MAIN  MENU 

IF(IMENU2  .EQ.  S)THEN 
Go  TO  5566 
ENDIF 

A****************************************************************** 


END  1  END  OF  MAIN  PROGRAM  'BODE* 

A 


******************  SUBROUTINE  SECTION  *********************** 

************************************************************9^*** 


**********************************************************M*******#*4r9^4r 
******************************************************************** 
******************************************************************** 
SUBROUTINE  REAOFREQ(GFREQ.GLEVEL,NTABLE) 

* 

*  Written  by  A.  L.  Helinskl  US  TACOM  AMSTA-RY 

It 

*  This  subroutine  reads  the  frequency  table  used  by  the  time  domain 

*  simulation  (ACSL)  to  create  the  Input  signal  to  excite  the  transfer 

*  function. 


INPUT;  None  (File  Table  Name) 

OUTPUT;  GFREQ(Index)  Table  Frequencies  (Hz)  (Frequency  of  each  s1n  wave) 


6LEVEL 

NTABLE 

REAL  GFREQ(50) 
CHARACTER*30  HLE 


Amplitude  of  each  sine  wave  (1  value  for  all) 
Number  of  Sin  Waves  (Frequencies) 


WRITEfB.lOO) 

table?',/, 

+  '  Example:  freqs.dat*) 

READ (5, 140) FI LEN 
140  FORMAT? A30) 

OPEN! 10 , FI LE»FILEN, F0RM= ' FORMATTED* , SHARED, 

*  STATUS='OLD',ERR=hO) 

READ(10,400)6LEVEL 
400  FORMAT(/,610.4) 


NTABLEal 
800  READ(10.600)6FREQ(NTABLE) 
600  FORMAT (010. 4) 

1F(GFREQ(NTABLE)  .LT.  0.)i 
NTA8LE>NTA8LE-fl 

SOTO  800 
RIT£75,1701 

170  fORMATC  ERMR  OPENING  FII 
300  CLOSE(lO) 

RETURN 

END 


(G10.4) 

EQ( NTABLE)  .LT.  0.)G0T0  300 


}R  OPENING  FILE') 


SUBROUTINE  GENERATE(SINPUT.6FREQ.6LEVEL, NTABLE, T) 

Written  by  A.  L.  Hellnski  US  TACOM  AMSTA-RY 

This  subroutine  creates  the  sweep  wave  (Sum  of  sine  waves) 
for  creating  the  Input  signal  to  a  frequency  response 
analysis  to  determine  transfer  function  characteristics. 

INPUT:  6FREQ  Table  of  Frequencies 
GLEVEL  Amplitude  for  all 
NTABLE  Number  of  Sine  waves  (Frequencies) 

OUTPUT:  SINPUT  Sum  of  sine  wave  signal  (Signal  Generated) 
T  Corresondlng  Time  (Seconds) 


DIMENSION  GFREQ(60) 

SINPUTsfl. 

Pts3.141S92654 
PO  100  I !»1, NTABLE 

5INPUf=SlNPUT  +  GLEVEL*SIN(2*PI*T*GFREQ(II)) 

100  CONTINUE 
RETURN 

**e*******4t*^*********4t4rit*ik*4r**4t****W*4t*1lr****************4t*4t***W* 

******4.********^********#*i(.****W*************************W*4t**4r**4t 

A***************************************************************** 

SUBROUTINE  INTEP(ICHAN.NSAMP,OELTA_F. FREQ, DB, PHASE, DATA) 

*  Written  by  A.  L.  Hellnski  US  TACOM  AMSTA-RY 

* 

*********************  INTERPOLOLATION  ************************ 

*  This  subroutine  Interpolates  points  from  the  frequency  response  results 
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* 

* 

* 

* 

* 

* 

* 


iGaIn  and  Phase)  to  re-astabllsh  data  with  a  constant  step 
Delta  Freq.  Hzj  so  that  the  further  processes  will  be  adaptible. 
n  other  words  this  prograw  basically  converts  a  set  of  raw  data 
of  varied  steps  (Varied  Saapling  Rates)  to  a  set  of  data  with  constant 
step.  If  you  think  that  sounds  hairy,  wait  until  you  see  the  data! 
Warning  must  be  a  sufficient  nuwber  of  points  to  begin  with,  enter 
at  your  own  risk.  The  wore  drastic  the  changes  In  the  data  the 
awre  points  will  be  required. 
*************************************************************************** 
* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


INPUT;  ICHAN 

NSANP 
DELTA  F 
FREO  “ 

DB 

PHASE 


Index  for  channel  ICHAN  will  be  gain 

and  ICHAN'4’1  will  be  phase 

Number  of  samples 
Desired  Step  frequency  (Hz) 

Actual  frequency  of  points  for  Phase  G  DB 
before  Interpolation.  (Un- Interpolated  points) 
dB  values  of  un- Interpolated  points 
Phase  (Deg)  values  of  un-interpolated  points 


OUTPUT:  NSAMP 

DATA(ISAMP,ICH) 


New  number  of  samples  for  Interpolated 
results  (Changed  from  Input) 

New  Interpolated  data  with  constant  step 


Generate  slopes  and  Intercepts  from  the  frequency  selected 
points. 

OIHENSION  PHASE(SO},0B(5O).FREQ(5O) 

OIHENSION  DB  SL(50).DB  ITISO) 

DIHENSION  PHXSE  SL(S0).PHAsE.IT(50) 

REALM  DATA!  10000,6) 

INTEGER*4  NSAMP 

00  1999  J«1,NSAMP-1 
PHASE_SL(J)  * 


(PHAOEI 


PHASE_It(J)  « 

ij-teEj 


)  -PHASE( J) )/(FREQ( J+1 ) -FREQ( J) ) 
SL(J)*FREQ(J) 


1999 


* 

* 

* 


PHASf(J 
OB  SL(J) 

(1JB(J+1)-0B(J))/(FREQ(J+1)-FREQ(J)) 

OB  lT(J)  * 

0B(J)-6b  SL(J)*FREQ(J) 

CONTINUE 

FREQFIN=FREQ(NSAMP) 

DETERMINE  NEW  INTERPLOATED  DATA 
NSAMP=INT(FREQFIN/0ELTA  F-1) 

step=deltA_f 

J=1 

DO  7510  Ksl, NSAMP 
FREQN=DELTA  F*(K-1) 

FREQ0=FREQP) 

IF(FREQN  .GT.  FREQO)THEN 
J=J+1 
ENDIF 

DATA(K,ICHAN)-0B  SL(0)*FREQN-«^DB  IT(J) 

0ATA(K,ICHAN+1)=PHASE  SL( J)*FRE0N+PHASE  IT(J) 

CONTINUE 
RETURN 
END 

************************************************************************** 


7510 
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*  «  «  « 


************************************************************************** 

**************************************************************************** 


SUBROUTINE  TEK  DELAY 
*  “ 

*  Written  by  AL  Reid  US  TACCM  ANSTA-RV  for  plotting  routine 


*  This  subroutine  w111  delay  2  seconds  to  allow  the  tektronlx  screen  ample 

*  time  to  clear  Itself. 

* 

*  INPUT:  NONE 

*  OUTPUT;  NONE 

* 

TNOW  «  SECNDS(O.O) 

10  DELTA  a  SECNOS(TNOW) 

IF  (DELTA  .LT.  2.0)  GOTO  10 

RETURN 

END 

* 


* 

SUBROUTINE  CROSS  OVER! ICH.NSAMP, DELTA  F, DATA) 

*****MM*«r**  determine  all  crossover  points  ♦*****^*^**^»»** 
* 

*  Written  by  A.  L.  Hcl Inski  US  TACON  ANSTA-RV 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


* 


This  subroutine  determines  the  crossovers  of  a  frequency 
response  data  for  stability  checks  for  the  ease  when  the  response 
Is  a  total  open  loop.  The  process  simply  scans  the  DATA  array 
for  magnitude  (DB)  and  Phase  and  detects  any  crossover  points  of 
0  dB  or  */-  IBO  Deg  respectively. 

INPUT;  ICH  Index  channel  as 

ICH  would  be  Gain 
then  ICHl  would  be  Phase 

NSAMP  Number  of  samples 

DELTA  F  Step  frequency  in  Hz 

0ATA(^AMP, Channel  Number)  Data  as  Gain  or  Phase 

OUTPUT;  NONE  (Results  are  printed  on  screen) 

INTERNAL;  SIGN  0B,SIGN_PH  Change  of  sign  Indicators 

JCRP.JCRG  Index  for  crossovers 

1NTE6ER*4  NSAMP 

DIMENSION  SIGN  DB( 10000). SIGN  PHI 10000) 

DIMENSION  JCRPT10000),OCR6(10DOO) 

REALM  OATA(10D00.6) 


*  DETERMINE  ALL  POINTS  AROUND  0  dB 

* 


DO  8199  J-1, NSAMP 


DETERMINE  ALL  POINTS  AROUND  0  dB 
IFjOATA^J^ICH)  .EQ.  0.)THEN 
ELSE  " 
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SIGN  OB(J)*OATA(J.ICH)/ABS(DATA(J,ICH)) 
EN01F“ 


DETERMINE  ALL  POINTS  AROUND  -180  DEG 

IF{DATA(J,ICH-»-l)  .EQ.  -180.)THEN 
SIGN  PH(J)=0. 

ELSE  “ 

IF(DATA(J,ICH+1)  .LT.  O.ITHEN 


IF(DATA(J,ICH+1)  .LT.  O.ITHEN 
SIGN  PHh)s(-i80.-DATA{j.ICH+l))/ABS(-180.-DATA(J,ICH-i-l)) 
ELSEIFTDATA{J,ICH+1)  .GT.  O.JTHEN 
SIGN  PH(J)«(-180.+DATA(J.lCHt-l))/ABS(-180.-t>DATA(J.ICH-i-l)) 


ENDIF 

ENDIF 

* 

8199  CONTINUE 
* 

*  DETERMINE  SIGN  CHANGES 

it 

ICG=1. 

ICP=1. 

DO  8198  J«l,NSAttP-l 


uw  n^nnr  A 

IFISIGN  08(J^1>  .ME.  SIGN  DB(J)>THEN 
JCR&IlC6i»j4-l 
IC6«lC&»i 
ENDIF 

IF(SIGN  PH(J>1)  .NE.  SIGN  PH(0))THEN 
JCRTfflCP)=J+l 
ICP«lCP+l 
ENDIF 
CONTINUE 

WRITE (5, 2005) 

FORMAT! MX.'********  GAIN  CROSSOVERS  ******** 

lOX,'  GAIN(dB)',13X,'  PHASEIDeg) ' , 13X, *  FREO(Hz)') 
I F ( I CG  . EQ.  1)THEN 
MRITE(5,9373) 

FORMAT ('**  NONE  **') 

ELSE 

DO  1666  J«1,IC6-1 


IV  AVVV 

FHZH=(JCRG(J))*DELTA  f 
FHZL=(JCRG(J)-iy*OELTA  F 
HRITE(5,4502jOATA{JCRGTJ)-l,ICH),DATA(JCR6(J)-l.ICH+ll 

.FHZL 

WRITE(S,4502}DATA(  JCRG!  J) ,  ICH)  .DATA!  JCRG(  J) .  ICH-t-1)  ,FHZH 
FORMAT(6X,Fl6.4,m,F10.4,13X,F10.4) 

WRITE (5,4903} 

FORMAT (/) 


FORMAT!/) 

CONTINUE 

ENDIF 

WRITE!5,447) 

READ(5.*) 

WRITE (5. 16671 
FORMAT !//,23X,' 


format !//,23X, *********  PHASE  ♦/-  180DE6  CROSSOVERS 
/.10X,»  PHASEjOE6)M3X,*  GAIN(DB) ’ ,13X,  •  FREQ(Hz)’) 
IF(ICP  .EQ.  1)THEN 


*******  • . 
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NRITE(S,9365) 

FORMAT(/,'  **  NONE  **  *) 

ELSE 

DO  1234  J=1.ICP-1 
FHZH=(JCRP{J))*OELTA  F 
FHZL=(JCRP(J)-mDELTA  F 

WRrTE(5.45d2)DAfA(JCRPTJ)-l,ICH+l),DATA(JCRP(J)-l,ICHJ^ 

WRITE j5, 4502 jOATA(JCRP(J),ICH+l),DATA(JCRP(J)iICH)lFHZH 

CONTINUE 
ENOIF 


HIT  RETURN  TO  CONTINUE*) 


return 

end 

*^*i*****l^*****ik**«******************************ik***ii^********ik****4r** 

^  SURROUTINE  ENTER_NPOLYS(N_MAT,NOOR,NPLyS.N_DRDER) 

*  Written  by  A.  L.  Hellnski  US  TACOH  AHSTA-RV 

*  this  subroutine  will  let  you  enter  the  Mlynbeileil. representing 

*  the  nueeretor  of  e  trensfer  function.  It’s  no  olffiTint  then 

*  Subroutine  ENTER^DPOLYS. 

*  INPUT;  NONE  (Enter  ell  output  penneeteirl  by  hllriU) 

*  OUPUT:  N  MAT (Poly#, ORDER#)  Polynd*ilBl  Cbefflclfents 

*  NHORiPOLY#)  Order  bf  each  poiynoeilel 

*  NPLYS  Number  of  PolynbiUI 

*  N  Order  Total  Order  (Degree  of  product) 

*  "  Sum  bf  all  DODR(POLYl) 


REAL  N  MAT(20.20) 
INTEGER  NOOR( 20). NPLYS 

WR1TE(5,10) 

F0RMAt(//////////////// 


I  NUM.  MATRIX  COEFlClENtS(POLY#,ORDER) 
I  NUM.  ORDER  (POLY#)  «  NuM.  #  OF  POLYS 


FORMAT  (//////////////////////////////////////////// , 

'  Enter  number  of  polynomials  In  the  numerator?') 
READ(5.*)NPLYS 
ilRITE(S.20) 

FORMAT (//) 

N  ORDER=0 


N  ORDER=0 
DO  90  1=1, NPLYS 
HRITE(5,30)1 

E0RMAt('  Enter  order  of  poly  #*,12) 

EAO(5,*)NOOR(I)  ^ 

N  OROER=N  OROER-fNOOR(I)  f  ADO  TOTAL  ORDERS  (DEGREE) 
WRITE  (5, 40) 

FORMAT (/) 

00  50  J=N00R(I)-M,1,-1 
WRITE(5,60)J-1 

FORMATC  S***,I2.*  TERM  IS  ?  *,:) 

REA0(5,*)  N_MAT(I,J)  I  NUMERATOR  MATRIX  (POLY#, ORE 
CONTINUE 
WRITE(5,70) 

FORMATC/) 

CONTINUE 


ilX  (P0LY#.0R0ER4-1) 
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RETURN 

END 

************************************************************** 
************************************************************** 
^  SUBROUTINE  ENTER_DPOLYS(0_MAT,DODR,DPLYS,D_OROER) 

*  Written  by  A.  L.  Hel Inski  US  TACOH  AMSTA-RY 


This  subroutine  will  let  you  enter  the 
the  denonlnator  of  e  transfer  function. 
Subroutine  ENTER  NPOLYS. 


polynoMlals  representing 
It's  no  different  then 


INPUT: 

OUPUT: 


NONE  (Enter  all  output  paraneters  by  hand) 
DMAT(POLY«. ORDER!)  Polynoalal  Coefficients 
OODRfPOLY!)  Order  of  each  polynomial 

OPLYS  Number  of  Polynomial 

D  Order  Total  Order  (Degree  of  product) 

Sum  of  all  DdDR(POLYI) 


dents 

ynomlal 


REAL  D  HAT (20 
INTEGER  D00R( 


.DPLYS.D_ORDER 


WRITE! 5» 10) 

FORMAtC/////////////////////////////////////, 

'  Enter  number  of  polynomials  In  the  denominator?') 
READ(5.MDPLYS 
WRITERS, 20) 

FORMAT!//) 

0  0RDER=0 
00  90  Isl. OPLYS 
WRITE(S,30)I 

F0RMAt('  Enter  order  of  poly  !  *,I2) 

READ(5.*)000R(I) 

0  ORDERED  ORDER-i^OOORd) 
write (5, 40) 


write  (5, 4(1 
FORMAT!/) 
00  50  J»0G 


00  50  J»OOOR{I)-t-l,l.-l 
WRITE(5,60)J-1 

60  FORMATC  S**',I2,'  TERM  IS  ?',:) 

REA0(5,*)  0  HAT(I.J)  I  NUMERATOR  MATRIX  (POLVf.ORDER-i-l) 
50  CONTINUE 

MRITE(5.70) 

70  FORMAT!/) 

90  CONTINUE 

* 

RETURN 

END 

******************************************************************** 
******************************************************************** 
SUBROUTINE  MULTIPLY  POLYS 

+  {T0TAL_0RDER,R0LY_0RDER,  NPOLYS, MATRIX,  PRODUCT) 

*  Written  by  A.  L.  Hellnski  US  TACOM  AMSTA-RY 

*  This  subroutine  multiplies  polynomials  for  resulting  product  of 

*  a  single  resulting  polynomial.  The  procedure  Is  much  the  same  way 

*  as  one  would  do  by  hand.  (Term  by  term) 

*  INPUT; 

* 


TOTAL  ORDER 


Total  degree  of  product  was  done  before  hand 
by  adding  the  degree  of  each  polynomial 
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POLY  ORDER  Degree 
NPOLTS  Nunber 
MATrjX( POLY#, Order#)  Matrix 


*  OUTPUT: 

* 


PRODUCT (Order  #) 

INTEGER  TOTAL  ORDER 
INTEGER  POLY  aRDER(20} 
INTEGER  NPOLYS 
REAL  HATRIX(20.20) 

REAL  SN(SO) 

REAL  PRd0UCT(40) 


Degree  of  each  Individual  po1yno«1a1 
Nuaber  polynomials 
Matrix  of  Polynoalal  coefficients 
(Slaply  a  2-D  array) 

Final  product  polynomial 

I  ORDER  OF  PRODUCT  (OR  TOTAL  SUN  degrees) 
!  ORDER  OF  EACH  POLY 
I  NUNBER  OF  POLYS 
t  COEFFICIENTS  OF  (POLY#,OROER} 

I  USED  AS  SUMNER  DURING  OPERATION 
I  PRODUCT  RESULTS 


PRODUCT  must  have  a  Intitlal  value. 
Initial  product  will  be  1. 


In  this  program  the 


PROOUCT(l)=l. 


INITIAL  PRODUCT  WILL  BE  1. 


LEVL=0 

DO  10  N=l, TOTAL  ORDER^-l 

SM(N)sO.  ~l  SUMMER  USED  IN  MULIPLICATION 

10  CONTINUE 

* 

DO  20  Isl. NPOLYS 
00  30  f&l.LEVL+l 

DO  40  L*l,POLY  OROER(I)-i-l 
SM7K+L-i)=SMrK+L-l)  ♦  MATRIX(I,L)  *  PRODUCT(K) 

40  CONTINUE 

30  CONTINUE 

DO  so  Mai, TOTAL  ORDER+1 
PRODUCT (M)aSMTM) 

SMMaO. 

50  CONTINUE 

LEVLaLEVL-^POLY  OROER(I) 

20  CONTINUE 

* 

RETURN 

END 

********************************************************************** 
************«'*************************«Hkr****|lr«r***  ********************* 

SUBROUTINE  POLY  RESPONSE (N  DEG, POLY  N.D  DEG, POLY  D, 

+  F  STEP,F  END,NrHAN  MAG,NCHAN_PHASE70A'n(,NPTS)  “ 

*  “  ~  ” 

*  Written  by  A.  L,  Hellnski  US  TACOM  AMSTA-RY 

*  Determine  frequency  response  for  a  single  polynomial  describing  the 

*  numerator  and  one  describing  the  denominator  of  a  transfer  function  for 

*  a  given  frequency  range. 


INPUT;  N  DEG 

POLY  N(Order#) 
D  DEO 

POLY  D{ Order#) 
F_STTP 

F  END 
NOHAN  MAG 
NCHAN“PHASE 


Degree  of  numerator 

Coefficients  describing  numerator  poly 

Degree  of  Denominator 

Coefficients  describing  denominator  poly 

Desired  delta  frequency  of  range 

Also  first  frequency  point  to  evaluate 

End  frequency  of  range  (Span) 

Desired  channel  number  for  Magnitude  (dB) 


Magnitude 


Desired  channel  number  for  Phase  (Deg 


* » *  »  »  »  » 


* 

* 

* 

* 

* 

* 

* 


* 


* 


OUTPUT;  NPTS  Number  of  frequency  pts  made 

OATA(NSAHP.NCHAN)  DATA  made 

where  Gain  (dB)  Is  DATAfNSAMP.NCHAN  NAG) 
Phase  (Deg)  Is  DATA(NSANP,NCHAN~PHASE) 
B  frequencies  «  NSAMP*F  STEP  Hz~ 


REAL  OREL, DING, NREL.NING 
REAL  PHASEN, PHASED 
REALM  data!  10000, 6) 

REAL  POLY  0(40), POLY  N(40) 

INTEGER  0_0EG.N_DE6 

PIM.141S926S4 

SCALsl80./PI 

EPSILON»10.E-30 

NPTS* I NT ( F_END/F_STEP ) 

DETERMINE  GAIN  AND  PHASE  FOR  NUMERATOR  AND  DENOMINATOR 

DO  90  JJ*1.NPTS 
FR*F  STEP*JJ 
WR*FR*2.*PI 

Determine  Real  and  Imaginary  components  for  NUN  and  DEN 
CALL  POLY  REAL  IMAGCNREL.NIHG.WR.N  DEG, POLY  N) 
CALL  P0LYI!REALIIMAG(DREL.DIH6,WR,DIDEG.P0LYID) 


Determine  dB  Gain 


* 


* 

* 

«t 

* 

* 


* 


GAIN»SQRT(  (NREL**2+NIHG**2)  /  (DREL**2+DIMG**2)  ) 

IFIGAIN  .LE.  EPSILON}THEN 
OATA(<]J,NCHAN  MAG)  *  -580. 

ELSE 

0ATA(JJ,NCHAN  MAG)  =  20.*  AL0610(GAIN) 

ENOIF 

Determine  PHASE  for  NUN  and  DEN 
Numerator 


IFfNREL  .EQ.  0.  .AND.  NIMG  .LT.  0.)THEN 
PHASEN  =  -90. 

ELSEIF(NREL  .EQ.  0.  .AND.  NIMG  .GT.  0.)THEN 
PHASEN  *  90 

ELSEIF(NIMG  .EQ.  0.  .AND.  NREL  .GT.  0.)THEN 
PHASEN  =  0. 

ELSEIFFnIMG  .EQ.  0.  .AND.  NREL  .LT.  0.)THEN 
PHASEN  *  180. 

ELSE1F(NREL  .GT.  0.  .AND.  NIMG  .GT.  0.)THEN 
RATION=ABS(NIMG/NREL) 

PHASEN*  SCAL*ATAN(RATION) 

ELSEIFInREL  .LT.  0.  .AND.  NIMG  .GT.  0.)THEN 
RATION=ABS(NIMG/NREL) 

PHASEN=180.  -  SCAL*AfAN(RATION) 

ELSEIFINREL  .LT.  0.  .AND.  NIMG  .LT.  0.)THEN 
RATI0N=A8S ( N IMG/NREL ) 

PHASEN=180.  +  SCAL*AtAN(RATION) 

ELSEIF(NREL  .GT.  0.  .AND.  NIMG  .LT.  0.)THEN 


RAT ION»ABS ( NIHG/NREL) 
PHASEN>360.  -  SCAL*AfAN(RATION} 


ENDIF 

DENOMINATOR  PHASE 


IFfDREL  .EQ.  0.  .AND.  0IM6  .LT.  0.)THEN 
PHASED  ::  -90. 

ELSEIFIDREL  .EQ.  0.  .AND.  DIHG  .6T.  0.)THEN 
PHASED  »  90. 

ELSEIFIDIHG  .EQ.  0.  .AND.  DREL  .6T.  0.)THEN 

PHASb  >  0. 

ELSEIFIDIHG  .EQ.  0.  .AND.  DREL  .LT.  0.)THEN 
PHASED  «  180. 

ELSEIFIDREL  .GT.  0.  .AND.  DIHG  .6T.  0.)THEN 
RATldD=ABS(DIHG/DREL) 

PHASED^  SCAL*ATAN(RATIOD} 

ELSEIFIDREL  .LT.  0.  .AND.  DIHG  .GT.  0.)THEN 
RATI0D=ABS(DIMG/DREL) 

PHASED^IBO.  -  SCAL*AfAN(RATIOD) 

ELSEIFIDREL  .LT.  0.  .AND.  DING  .LT.  0.)THEN 
RAT lOO-ABS ( DING/DREL) 

PHASED>180.  -t-  SCAL*ATAN(RAT10D) 

ELSEIFIDREL  .GT.  0.  .AND.  DING  .LT.  0.)THEN 
RATI0D=ABS(DIHG/DREL) 

PHASED«360.  -  SCAL*ATAN(RATIOD) 


ELSEIFII 

RATIOI 


0.)THEN 


0.)THEN 


ENDIF 

Determine  Totel  Phe$e 


DATA ( J J , NCHAN.PHASE )=PHASEN- PHASED 

90  CONTINUE 

* 

RETURN 

END 

*********************************************************************** 

*********************************************************************** 

*********************************************************************** 

^  SUBROUTINE  POLY_REAL_IMAG(XREAL,XIMAG,OMEG,NORDER,POLY) 

*  Written  by  A.  L.  He! Inski  US  TACOM  AMSTA-RY 

* 

*  DETERMINE  REAL  AND  IMAGINARY  COMPONENTS  FOR  A 

*  GIVEN  FREQUENCY (OMEG)  AND  POLYNOMIAL 

*  This  subroutine  called  upon  by  sub-poly-response  to  determine  the 

*  Real  and  Imaginary  components  for  a  given  polynomial  and  frequency(W). 

*  It  simply  plugs  In  S=JW  and  resolves  the  complex  algebra.  (  W  omega  In 

*  rad/sec  and  J  Is  complex  where  j^SQRT  [  -1  ].  ) 

* 

*  INPUT:  POLY  (Array  of  POLY  coefficients) 

*  NORDER  Order  of  polynomial 

*  OMEG  Frequency  In  Rad/Sec 

*  OUPUT;  XREAL  Real  Component 

*  XIMAG  Imaginary  Component 
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REAL  POLY (40) 

XREAL-0. 

XIHA6=:0. 

DO  10  II=1.N0R0ER4^1 

I0R0=II-1 

INDEX:>I0RD 
IF( INDEX  .GE.  4) 
INDEX»INDEX-4 
60  TO  20 
ELSE 


IF(1N0EX  .EQ.  0)THEN 
XREAL>XREAL  *  POLYflORO^-l) 
ELSEIFIINDEX  .EQ.  ijTHEN 
XINAG=XIMA6  *  POLY(IORD^l) 


*0ME6**I0RD 


XIMAG  ♦  POLY(IORD+1)*OMEG**10RD 

NOEX  .EQ.  2) THEN 

XREAL  -  P0LY(I0RD^1)*0MEG**I0RD 

NDEX  .EQ.  3) THEN 

XIMAG  -  P0LY(I0RD+l)*0ME6**10RD 


ELSEIF(INOEX  .EQ.  2) 
XREAL<XREAL  -  POLY 
ELSEIFFinDEX  .EQ.  3) 
XIHAG«XIHA6  -  POLY 
ENDIF 
ENDIF 


10  CONTINUE 

* 

RETURN 

END 

********************************************************************* 
********************************************************************* 
********************************************************************* 
SUBROUTINE  LIST  POLY(IO,POLY  MATRIX, PRODUCT  POLY, 

^  ♦  0RDER_MAT11IX,NUM_P0LYS, DEGREE) 

*  Written  by  A.  L.  Hel Inski  US  TACOM  AMSTA-RY 

*  This  subroutine  list  the  polynomial  coefficients. 

*  INPUT;  ID  Identify  Individual  polnomlals  or  sinale  product 

*  If  ID  s  1  Number  of  Individual  Polynomials 

*  ID  =  Z  Single  Product  Polynomial 

*  POLY  MATRIX(POLY#, Order  #) 

*  Actually  describes  the  coefficints  In  this 

*  2-D  array  fashion.  (When  10=1  otherwise  a  dummy) 

*  PORDUCT_POLY(Order  #) 

*  Oecribes  a  single  polynomial  In  coefficient  terms. 

*  (When  ID=2  otherwise  a  dummy  variable.) 

*  OROER_MATRIX(POLY  #) 

*  Describes  the  degree  for  each  polynomial 

*  NUM_P0LYS  Number  of  polynomials  (=1  for  ID=2) 

*  DEGREE  Total  degree  of  polynomial  (or  sum  of  orders) 

*  OUTPUT:  None  (Just  simply  lists  the  polynomials  on  the  screen) 


10 


20 


30 


* 

* 

* 


70 


80 


* 

90 


REAL  MATRIX(20.20) 

INTEGER  lO.NUH  POLYS.OEG, ORDER (20) 
CHARACTER*!  RfPLY 


WRITEf  .  , 

FORMAt ( /////////////////////////////////////////// , 

*  '  Do  you  desire  to  re-enter  the  entire  NUMERATOR  ?  (Y  or  N)'] 

READ(5,20)REPLY 
FORMAT  A) 

IF(REPLY  .EQ.  'Y*  .OR.  REPLY  .EQ.  •y')THEN 

CALL  ENTER  NPOLYS(MATRIX, ORDER, NUM  POLYS.OEG) 

RETURN 

END  IF 
ELSE 

WRITERS. 30) 

FORMAT ( ////////////////////////////////////////////. 

*  '  Do  you  desire  to  re-enter  the  entire  DENOMINATOR  ?  (Y  or  N)') 
REA0(S,20)REPLY 

IF(REPLY  .EQ.  'Y*  .OR.  REPLY  .EQ.  •y')THEN 

CALL  ENTER  0P0LYS(MATRIX, ORDER, NUM  POLYS.OEG) 

RETURN 

ENOIF 

ENOIF 

Edit  Polynoalal  by  terns 

IF(NUM  POLYS  .GT.  1)THEN 
WRITE (5. 70) 

fOl>^kn////////////////////////////////////////, 

*  I  Wf*!?**  Polynomial  do  you  desire  to  change  ?I*^» 


reI 


Give  the  number  representing  the  polynomial  sequence  )') 
READ(5,*)MP0LY 


ELSE 

MP0LY*1 

ENOIF 

IF(0R0ER(MP0LY)  .GT.  1)THEN 
WRITE(5,80)MP0LY 

format I ///////////////////////////////////////// , 

+  '  Which  ORDER  coefficient  of  S  do  you  desire  to  •,/, 

♦  •  change  In  POLY  #  M2,'  ?’) 

READ(5,*)MC0EFF 

ELSE 

MC0EFF=1 

ENOIF 

WR ITE ( 5 , 90 ) MATR I X ( MPOLY , MCOEFF+1 ) 

formaT(//////////;///////////////////////////////////. 

+  •  Change  ',F10.4,'  TO  ?•) 

REA0( 5 , * ) MATR I X ( MPOLY , MCOEFF+1 ) 


RETURN 

END 

SUBROUTINE  READERO(TIME, SIGNAL. NSAMP, 
+  CHANNEL, UNITS, NCH) 
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* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


Written  by  AL  Reid  US  TACOH  AHSTA-RY 
Modified  by  A.  L.  Hel Inski  US  TACOM  AMSTA-RY 

This  subroutine  orignally  written  by  A1  Reid  and  Modified  by  A.  L. 
Hellnskl.  It  simply  reads  a  data  file  of  a  time  history  In  ERD  file 
format.  For  this  application  ACSL  creates  the  time  data  file 
simulating  the  output  of  a  transfer  function  this  routine  reads  the 
time  history  Into  the  desired  variables  by  calling  It. 


INPUT: 

OUTPUT; 


None  (Data  File  Name) 

SIGNAL  Amplitude  of  history 

Time  corresponding  Time  Is  Seconds 


DIMENSION  TIMEaOOOOl, SIGNAL!  10000) 

CHARACTER*80  ERD  TITLE. LONG, TITLE, DUMHYSO 
CHARACTER*64  ERO  FILE.HOR  FILE, ERD  FILE  O.HDR  FILE  0 
CHARACTER*32  LONB  NAME ( 6 )70UMMY32, CHANNEL 
CHARACTER*12  DUMMY 

CHARACTER'S  SHORT  NAME(6).UNIT  NAME(6).XUNIT,0UMMY8. 

+  UNITS 

CHARACTERM  ERD.HDR 

CHARACTER*!  COMMA.REPLY 

REAL*4  SCALE(6] .0FFSET(6) .ROATA(6. 10000) 

INTE6ER*4  NSAMP 

INTE6ER*4  START  SAMP  ELIN.END  SAMP  ELIM 
INTEGER*2  ERO  URIT.HDR  UNIT 
DIMENSION  IAR11AY(400),TDATA(4) 

DATA  ERO  UNIT.HOR  UNIT/10, 11/ 

ERO  «  '.FRO' 

HOR  =  ’.HOR’ 

* 

*  Determine  the  name  of  the  Input  data  file 
10  HRITE(5,20) 

20  FORMAT!////,*  Enter  name  of  ERO  data  file  to  analyze?') 

READ! 5. 30)  ERD_FILE 


format{a; 
CALL  str: 
HDR  FILEi 


TRIMIHDR  FILE.ERO  FILE.LENGTH) 
LENGtH+lTLENGTH+4T  *  HOR 


ERD"FlLE(LENGTH+l;LENGTH+4)  »  ERD 

*  ~ 

*  Open  the  data  file,  print  the  header  characteristics,  and  determine  If 

*  this  Is  the  correct  data  file 

* 

0PEN!HDR  UNIT,FILE=HDR  FILE,FORM=*FORMATTED' , 

^  ♦  5HARE0,STATUS=''"0L0',ERR=210) 

*  Read  the  header  data 

* 

REA0!H0R  UNIT, 60)  DUMMY 
60  F0RMATIAT21 

READIHOR  UNIT, 70)  ERD  TITLE 
70  FORMAT (A80) 

REA0!H0R  UNIT.80)  NCHAN, COMMA, NSAMP, COMMA, NLINES.COMMA.NBIN, 

+  COMMA  .NBYTE ,  COMMA,  KEYNUM,  COMMA ,  STEP ,  COMMA ,  KEYOPT 

80  FORMAT!6!I7,A),E13.6,A,i7) 


READ{HOR“UNIT 

REA0!H0R”UNIT 


.100)  ! 
,110  ( 
.100  ( 


LONG  NAME(L).L:=1, NCHAN) 
UN I TINAME ( L ) , L» 1 , NCHAN ) 
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F0RHAT(18(E13.6,A)) 
FQRMAT(ai|A8|j 


90 

no  F0RNAT(7(A32] 

*  Write  out  the  header  Information 

* 

WRITERS, 120)  ERO  TITLE. NCHAN. NS AHP, STEP 
120  FORMAT(//,*  The  fitle  for  this  file  Is:*,/,*  *,A80,//,*  There  are 
■•>*,12,*  channels  of  data.',/,*  There  are  ',17,*  samples  for  each  da 
■fta  channel.*,/,*  The  step  size  Is  ',F8.5,'  seconds.',//) 

*  See  If  there  are  more  than  16  channels  to  plot 

IF  (NCHAN  .GT.  16)  THEN 

Type*,*  There  are  more  than  16  channels  to  read.* 

TYPE*,*  Please  FORGET  ABOUT  IT* 

CLOSE  (HDR  UNIT) 

STOP 

ENDIF 

* 

*  Write  out  the  additional  descriptor  lines 
IF  (NLINES  .GT.  0)  THEN 

Type*. 'The  following  are  the  optional  descriptor  lines:' 

DO  130  L>1.NLINES 

READCHDR  UNIT. 70)  LONG 
WRITE(S.T25)  long 
FORMATC  *,A80) 

CONTINUE 


125 

130 


ENOIF 


*  Is  this  the  correct  data  file 

it 

HRITE(5,140) 

140  FORHATu/, *ils  this  the  correct  data  file  to  analyze  (y  or  n)? 

*r{aO(5,150)  reply 
formatIa) 

IF  (REM  .EQ.  *N*  .or.  reply  .EQ.  *n*)  THEN 
CLOSE(H0R  UNIT) 

WRITEiS.lBO) 

vou  wish  to  look  at  another  file  (y  or  n)?  *) 


150 


160 


WRITE (5, 160) 
FORMAT (/.*$6o  you 
READ(5.1S0)  RE^Y 


IF  (REPLY  .EQ. 
GOTO  10 
ENDIF 


*N*  .OR.  REPLY  .EQ.  *n*)  STOP 


*  Open  the  data  part  of  the  file 

* 

IF  (KEYNUM  ,EQ.  5)  THEN 

0PEN(ERD  UNIT,nLE=ERD  FILE, STATUS* 'OLD* 
*  , SHARED, FORM*  *  FORMATTED* ) 

ELSE 

0PEN(ERD  UNIT,FILE=ERD  FILE,STATUS**OLD* 
,SHARED,FORM=*UNFORMATTED*) 

ENOIF 

CLOSE(HOR_UNIT) 

*  Read  the  data 

* 

J  *  0 
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170 

IF  (KEYNUM  .EQ.  5)  THEN 


REAOIERD  UNIT, 180. ERR=>220.END»230)  (R0ATA(I,J).I«1.NCHAN) 
F0RHAT(19(E13.6)) 


IF  (KEYNUM  .EQ.  0)  THEN 

READIERO  UNIT.ERR>220,EN0=>230)  (I0ATA(I).I«1,NCHAN) 
DO  idO  K^l.NCHAN 

RDATAJK^J)  »  FLOATJ(IDATA(K)) 

ELSE 

READ(ERO_UNIT. ERR=220 .END>230)  (RDATA(I.J),I«1.NCHAN) 
END  IF 


ENOIF 
GOTO  170 


210  TYPE*, ’Error  opening  data  file* 

STOP 

* 

220  TYPE*. ’Error  reading  data  In  file’ 

230  close! ERD.UNIT) 

* 

***********  DATA  IS  READ  IN,  NOW  START  EVALUATION  ******* 

* 

***  Convert  UNscaled  and  UNblased  data  to  proper  values 
DO  6002  I=1.NCHAN 
DO  6003  J=1.NSAMP 

RDATA(I,J)  =  RDATAd.J)  /  SCALE(I)  +  OFFSET(I) 
6003  CONTINUE 

SCALE(I)=1. 

OFFSET(1)=0. 

6002  CONTINUE 


NRITE(5,1018} 

FORMAT!//, ******  CHANNELS  *******',//) 

DO  1050  JI=1,NCHAN 

IF(JI  .EQ.  6  .OR.  JI  .EQ.  12  .OR.  JI  .EQ.  18  .OR. 
JI  .EQ.  24  .OR.  JI  .EQ.  30  .OR.  JI  .EQ.  36)  THEN 
NRITE(S,1032) 

READ{5,*) 

ENDIF 

WRITERS. 1060} JI, LONG  NAME(JI),UNIT  NAME(JI) 
FORMAT(/.’  CHANNEL  ’712,/, 1X,A32,/71X,A8) 

CONTINUE 

NRITE(5,1032) 

FORMAT (/, **********  HIT  RETURN  ********•) 
READ(5,*) 

WRITE(5,1081) 

FORMATC*  ENTER  CHANNEL  NUMOER  TO  ANALYZE') 
REA0(S,*)NCHANPS0 

UN I T5=UN I T_NAME ( NCHANPSD) 
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CHANNEL»LONG_NAME ( NCHANPSD ) 


* 

DETERMINE 

* 

DO  1061 

TIME(I 

SIGNAL 

1061 

CONTINUI 

* 

RETURN 

END 

Ja(N^I 


ROAtA(NCHANPSD,I) 


***********************************************************4'***** 

SUBROUTINE  TF  PLOT(ITERH.ERD  TITLE. NSAHP. STEP. NCHAN, DATA. 

+  L0N6_NAME7UNIT_NAME) 


Must  Insert  a  subroutine  for  plotting  ERD  data  format. 

This  particular  application  uses  a  PLOTIO  Library  routine 
(Not  Included).  It  Is  advised  to  use  a  subroutine  which  has 
a  graphics  mode  capable  of  plotting  semi-log  graphs. 

The  following  format  of  Input  variables  was  used. 


INPUT:  ITERH  10  for  terminal 

ERD  TITLE  Title  for  plot  ('NOTE*  from  main) 

NSARP  Number  of  samples 

STEP  Delta  F  (Hz) 

NCHAN  Number  of  Channels  (6  for  this  program) 

OATA(Nsamp, Channel  Number)  Data  In  amplitude 
LONG  NAME (Channel  Number)  Name  for  channel 
UNIT~NAME( Channel  Number)  Unit  for  channel 
OUTPUT:  None  (A  plot  on  the  screen!  what  else  do  you  expect! !l) 


ITERH 
ERD  TITLE 
NSARP 
STEP 
NCHAN 


•  TF^PLOT  Subroutine... 

* 

RETURN 

END 

* 

******  FOURI ER  SECTION  **************************************^* 
****************************************************************** 
****************************************************************** 
SUBROUTINE  FOURIER(X.Y.NPT. REAL. AIHAG, FREQ, MM) 

* 

*  Must  Insert  a  fourler  analysis  subroutine. 

*  The  one  used  for  this  analysis  was  developed  by  IEEE  for 

*  digital  processing.  (Not  Included) 

*  It  determines  the  fourler  coefficients  for  a  time 

*  history.  The  essential  requirement  Is  that  the  time  history  be 

*  described  by  a  Power-0f-2-Number-0f-Po1nts. 


INPUT; 


OUPUT: 


Y  (Nsamp) 
X  {Nsamp ) 
NPT 


Signal  In  Amplitude 
Corresponding  Time  (Seconds) 
Correspondlno  Number  Of  points  of  signal 
(Program  will  help  you  choose  closest 


(Program  will  help  you  choose  close 
power  of  2) 

REAL  Real  Coefficient  (Cos  Term) 

AIMAG  Imaginary  Coefficient  (Sin  Term) 

FREQ  Corresponding  Frequency  (Hz) 

MM  Number  of  points  In  frequency  domain. 


FOURIER  Subroutine... 


RETURN 

END 
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APPENDIX  C 

VARIOUS  STATES  OF  A  STEP  RESPONSE 


C-1 


STEP  RESPONSE  1  IN 

POSITION  COMMAND 


STEP  RESPONSE  1  IN 


RATE  COMMAND 


BTi] 


StEP  RESPONSE  1  IN 


FLOW  IN  (for  down) 

8 

- 


Flow  out  (tor  down) 


a 

o 


o 

a 

g 

a 

T 

og 

O 

4 

1111a. 

O 

o 

"^O.t 

M  0 

20  0 

.40  0 

.60  0 

.90  1. 

APPENDIX  D 

VARIED  STEP  RESPONSE  RESULTS 


VARIED  STEP  RESPDNlSe 


POSITION  COMMAND  .05  IN 

o  SIMULATED  ACTUATOR  RESPONSE 


POSITION  COMMAND  0.1  IN 

^  SIMULATED  ACTUATOR  RESPONSE 
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POSITION  COMMAND  0.25  IN 


VARIED  STEP  RESPONSE 


POSITION  COMMAND  0.5  IN 


SIMULATED  ACTUATOR  RESPONSE 


SIMULATED  ACTUATOR  RESPONSE 


VARIED  STEP  RESPONSE 

POSITION  COMMAND  0.75  IN 

8 


SIMULATED  ACTUATOR  RESPONSE 


APPENDIX  E 
VARIED  MASS  RESULTS 
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